Computation of dynamic adsorption with adaptive integral, finite difference, and finite element methods

Journal of Colloid and Interface Science 258 (2003) 310–321 www.elsevier.com/locate/jcis

Computation of dynamic adsorption with adaptive integral, finite difference, and finite element methods Ying-Chih Liao, Elias I. Franses, and Osman A. Basaran ∗ School of Chemical Engineering, Purdue University, West Lafayette, IN 47907-1283, USA Received 13 February 2002; accepted 28 October 2002

Abstract Analysis of diffusion-controlled adsorption and surface tension in one-dimensional planar coordinates with a finite diffusion length and a nonlinear isotherm, such as the Langmuir or Frumkin isotherm, requires numerical solution of the governing equations. This paper presents three numerical methods for solving this problem. First, the often-used integral (I) method with the trapezoidal rule approximation is improved by implementing a technique for error estimation and choosing time-step sizes adaptively. Next, an improved finite difference (FD) method and a new finite element (FE) method are developed. Both methods incorporate (a) an algorithm for generating spatially stretched grids and (b) a predictor–corrector method with adaptive time integration. The analytical solution of the problem for a linear dynamic isotherm (Henry isotherm) is used to validate the numerical solutions. Solutions for the Langmuir and Frumkin isotherms obtained using the I, FD, and FE methods are compared with regard to accuracy and efficiency. The results show that to attain the same accuracy, the FE method is the most efficient of the three methods used.  2003 Elsevier Science (USA). All rights reserved. Keywords: Dynamic adsorption; Numerical solutions; Integral method; Finite difference method; Finite element method; Diffusion-controlled adsorption; Surface tension

1. Introduction Dynamic surface tension behavior at fluid/fluid interfaces is important to many applications involving foams [1], coatings [2], lung surfactants [3], and drop and jet breakup [4]. The relation between the adsorbed surface density of surfactant and the dynamic surface tension at the interface permits studying interfacial transport by dynamic surface tension measurements. For a nonionic surfactant below its critical micelle concentration (cmc), the dynamic surface tension of a newly created surface is governed by a twostep process [5]. First, the surfactant monomers in the bulk solution diffuse to the subsurface layer. Then, the surfactant is adsorbed from the subsurface layer to the surface layer to affect the surface tension. To model the two-step process in a spatially one-dimensional problem, the diffusion in the bulk solution is described by a transient diffusion equation coupled with a time-dependent surface density. The general solution for the surface density is an * Corresponding author.

E-mail address: [email protected] (O.A. Basaran).

integral [6] relating the surface density and the subsurface concentration. To complete the mathematical statement of the problem, the relationship between the surface density and the subsurface concentration, which depends on the adsorption/desorption mechanism, is needed. Often, a local equilibrium assumption or a so-called “diffusion controlled” mechanism applies [6–9]. The adsorbed surface density is related to the subsurface concentration by a dynamic isotherm, which usually has the same functional form as the equilibrium adsorption isotherm [10,11]. Once the dynamic adsorption isotherm is specified, the problem is well posed and can be solved explicitly for the surface density and the bulk concentration profile. The problem of diffusion-controlled adsorption with a linear adsorption isotherm was first formulated and solved analytically by Sutherland [7]. However, since most of the isotherms are nonlinear for high surface tension reduction, numerical solutions to the problem are needed. The oftenused numerical method is applying the trapezoidal rule approximation to the Ward and Tordai equation [12,13]. Another approach is based on the idea that over small increments of time, the relationship between the adsorbed surface density and the sublayer concentration is nearly linear [14].

0021-9797/03/$ – see front matter  2003 Elsevier Science (USA). All rights reserved. doi:10.1016/S0021-9797(02)00096-6

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

This approach then builds a solution by linearization of the dynamic isotherm that appears in the Ward and Tordai equation over a large number of successive small intervals in time. However, an analytical solution with the linearized isotherm is then required at each step. Thus, the solution procedure is more complicated compared to the situation in which the isotherm is linear and not easily formulated. Finite difference (FD) methods for this problem have also been used to directly solve the differential equations governing the dynamic adsorption problem without using the Ward and Tordai method. For diffusion-controlled adsorption in a spatial domain of infinite extent, Miller developed an FD algorithm based on the Crank–Nicholson method [8]. The time-step size was variable and increased by 1% in each step. The spatial grid was uniform but empirically adjusted in time. For dynamic adsorption within a finite spatial domain, Chang and Franses developed a fully implicit FD scheme using the control volume method [11]. This scheme was further extended to pulsating area conditions in spherical coordinates [15]. Although several numerical methods for solving the onedimensional diffusion-controlled adsorption problem are available, the accuracy and efficiency of various methods have heretofore not been determined. To date, the convergence of the trapezoidal rule approximation to the Ward and Tordai equation could only be ensured by decreasing the time-step size. In the published FD schemes, previous authors either used small time steps of uniform size or empirically adjusted them, but did not guarantee the accuracy of computations. Clearly, what are needed are methods that incorporate algorithms for systematically generating spatial grids for minimizing the number of grid points used, employ adaptive time-integration techniques, and provide precise estimates of computational accuracy (or error). The aim of this article is to investigate the accuracy and efficiency of key computational methods for analyzing onedimensional diffusion-controlled adsorption. First the integral (I) method, which entails solving a Ward and Tordai type integral equation when the extent of the spatial domain is arbitrary, is put on a firmer foundation. The time truncation error of the trapezoidal rule approximation for the I method is evaluated, which allows time-step sizes to be chosen adaptively to ensure solution accuracy and computational efficiency. Next, an improved FD scheme and a new finite element (FE) approximation are developed. In both methods, a stretched grid is used for spatial discretization, which makes it possible to obtain mesh-insensitive solutions with drastically fewer grid points. Moreover, time integration with both the FD and FE methods is carried out using a predictor–corrector method and time-step sizes are determined adaptively. Finally, solutions obtained with the I, FD, and FE methods for nonlinear isotherms are compared with respect to their accuracy and efficiency. The new algorithms should be of interest for several reasons to researchers who work on problems of diffusion and adsorption of surfaceactive substance at interfaces. First, many researchers who

311

work on such problems use primarily the I and FD methods but without evaluating computational errors in or optimizing the efficiency of their algorithms. Second, the FE method, while a popular technique for solving ordinary and partial differential equations, has not been used in past studies of such problems.

2. Key equations and numerical methods 2.1. Governing equations For dynamic adsorption of a nonionic surfactant below its cmc, molecules are transported through a diffusion layer of length l and adsorb onto the interface. For a stagnant medium, l would be practically infinite. However, because of mixing, convection in the bulk phase can increase the rate of adsorption by decreasing the effective diffusion layer thickness and maintaining a nearly constant bulk concentration at x = l, where x is the coordinate measured from the interface. A model with a finite diffusion layer thickness was developed by Mysels and Frisch [16] and Chang and Franses [11] and will be briefly reviewed here. The one-dimensional diffusion equation in the x-direction is ∂ 2C ∂C = D 2 , 0 < x < l, t > 0, (1) ∂t ∂x where t is time, C(x, t) is the bulk concentration, and D is the bulk diffusivity, which is assumed to be constant. The uniform surface density Γ (t) is related to the concentration gradient near the interface x = 0 via the well-known mass balance equation,
∂C

dΓ =D , t > 0. (2) dt ∂x
x=0 The second boundary condition is that the concentration beyond the diffusion layer is fixed and equal to the bulk concentration C0 , C(l, t) = C0 ,

t > 0.

(3)

If the mixture is placed in a container so that the distance from the wall to the interface lc is small compared to the effective diffusion boundary layer thickness l, one may have to use the Neumann boundary condition ∂C/∂x = 0 at x = lc . Then, the calculated dynamic surface tension will equilibrate slightly faster than in the case of Eq. (3). For lc > l, the two results will be about the same, and for lc  l, the two results will be identical. Here, Eq. (3) is used throughout the article. The initial condition for C is C(x, 0) = C0 ,

0
x
l.

(4)

The initial condition for Γ depends on how the interface is created. If a new interface is created so rapidly at time zero that the surface is nearly clean, then Γ (0) = 0.

(5)

312

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

Because in Eqs. (1)–(5) C(x, t) and Γ (t) are coupled, an additional relation between Γ (t) and C(0, t) is needed. The relation depends on the mechanism of adsorption, whether the adsorption is diffusion-controlled, kinetic-controlled, or “mixed kinetics,” that is, affected by both diffusion and adsorption/desorption. Here, we examine the first case. We assume that the surface density Γ (t) and subsurface layer C(0, t) are in local equilibrium; i.e., adsorption/desorption is much faster than diffusion. The dependence of the surface density Γ (t) on C(0, t) is given by the dynamic adsorption isotherm, which is normally taken to have the same functional form as the equilibrium adsorption isotherm. The most often-used equilibrium isotherm is the Langmuir equation, Γe = Γm,L

KL C0 , 1 + KL C0

(6)

where Γe is the equilibrium surface density, Γm,L is the maximum surface density for the Langmuir model, and KL is the adsorption equilibrium constant. Typical values for Γm,L range from 8 × 10−6 mol/m2 (≈ 20 Å2 /molecule) to 1 × 10−6 mol/m2 (≈ 160 Å2 /molecule); for KL , the values fall in a much wider range spanning several decades, and typical values of KL C0 for nonionic surfactants at the cmc are about 100–300 [10]. The dynamic Langmuir isotherm is KL C(0, t) , Γ (t) = Γm,L 1 + KL C(0, t)

(7)

where Γm,L and KL are the same parameters as in the equilibrium Langmuir equation. The surface equations of state for the equilibrium and the dynamic Langmuir isotherm are   Γe , γe − γ0 = Γm,L RT ln 1 − (8) Γm,L   Γ (t) , γ (t) − γ0 = Γm,L RT ln 1 − Γm,L

(9)

where γe and γ (t) are the equilibrium and dynamic surface tensions, respectively, γ0 is the surface tension of pure solvent, R is the universal gas constant, and T is the absolute temperature. Another widely used equilibrium isotherm is the Frumkin isotherm, Γe = Γm,F

KF C0 , + KF C0

eAXe

(10)

where A is the binary interaction parameter, Γm,F is the maximum surface density, Xe ≡ Γe /Γm,F is the equilibrium surface coverage, and KF is the adsorption equilibrium constant. The dynamic Frumkin isotherm is KF C(0, t) , Γ (t) = Γm,F AX(t ) e + KF C(0, t)

Frumkin isotherms are   AXe2 , γe − γ0 = Γm,F RT ln (1 − Xe ) − 2    AX(t)2 . γ (t) − γ0 = Γm,F RT ln 1 − X(t) − 2

(12) (13)

For convenience in further mathematical treatment, the concentration, surface density, distance, and time are nondimensionalized as follows: C , C0 Γ θ≡ , Γe x ξ≡ , l DC02 t τ≡ . Γe2 u≡

(14) (15) (16) (17)

In the last equation, the characteristic time is taken as L2 /D, where L ≡ Γe /C0 is the adsorption length or the length of depletion [17], or the thickness of the solution containing an equal number of moles as in the equilibrium surface layer. The resulting dimensionless diffusion layer thickness is Nl ≡

lC0 . Γe

(18)

The smaller the value of Nl , the thinner is the diffusion layer, and hence the higher is the adsorption rate. The original dimensional Eqs. (1)–(5), when cast into dimensionless form, become ∂u ∂ 2 u = 2 , 0 < ξ < 1, τ > 0, ∂τ ∂ξ
∂u

dθ = Nl , τ > 0, dτ ∂ξ
ξ =0

Nl2

(19) (20)

u(1, τ ) = 1,

τ > 0,

(21)

u(ξ, 0) = 1,

0
ξ
1,

(22)

θ (0) = 0.

(23)

To complete the nondimensionalization, the dynamic isotherms are also normalized. The dimensionless dynamic isotherms are listed in Table 1 as 

θ (τ ) = G θ (τ ), us (τ ) , (24) where G is the functional form of each isotherm, and us (τ ) ≡ u(0, τ ) is the dimensionless subsurface layer concentration. Because the dynamic isotherms are nonlinear, Eqs. (19)–(24) need to be solved numerically. 2.2. Integral (I) method

(11)

where X(t) ≡ Γ (t)/Γm,F is the dynamic surface coverage. The equations of state for the equilibrium and dynamic

Following the analysis of Ward and Tordai [6], the general solution of Eqs. (19)–(24) has been obtained using the Laplace transform method, and is given by the integral

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

313

Table 1 Summary of isotherms used Henry

Langmuir

Frumkin

Equilibrium isotherm

KH C0

KL C0 Γm,L 1+K L C0

Equilibrium surface pressure

KH C0



−Γm,L RT ln 1 − ΓΓe

 AX 2  −Γm,F RT ln (1 − Xe ) − 2 e

KH

KL , Γm,L

KF , Γm,F , A

Dimensionless dynamic isotherm G(us , θ )

us

(1+NcL )us 1+NcL us

NcF us Xe (eAXe θ +NcF us )

Dimensionless surface pressure Π (θ )

θ

1−

Dimensionless parameters

–

NcL ≡ KL C0

Parametersa

K C Γm,F AXe F 0 e +KF C0 m,L

ln[1+NcL (1−θ )] ln(1+NcL )

ln(1−Xe θ )−(A/2)Xe2 θ 2 ln(1−Xe )−(A/2)Xe2

NcF ≡ KF C0 , A

a X is calculated from the equilibrium Frumkin isotherm. e

equation (in dimensionless form) 2 θ (τ ) = √ π

ζ =τ

where p  2  F (τi ) + F (τi−1 ) √ √ τi − τi−1 , I1 ≡ √ 2 π

F (ζ ) d ζ

i=1

ζ =0

2 +√ π

ζ =τ

p−1 

us (ζ )F (τ − ζ ) d τ − ζ ,

(25)

ζ =0

where F (τ ) ≡ 1 + 2

∞ 

e−(Nl /τ )n 2

2

n=1

and ζ is a dummy variable. In the limiting case of Nl → ∞, or as l → ∞, F (τ ) approaches unity and Eq. (25) reduces to the well-known Ward and Tordai equation: 

2 τ +√ θ (τ ) = 2 π π

ζ =τ

us (ζ ) d τ − ζ .

(26)

ζ =0

us (τi )F (τp − τi ) + us (τi−1 )F (τp − τi−1 ) 2 I2 ≡ √ 2 π i=1



× τp − τi − τp − τi−1 ,  (τp α≡ , π and

−us (τp−1 )F ((τp ) (τp β≡ √ . π Inserting the expression for θ (τp ) in Eq. (28) into the dimensionless dynamic isotherms G(us , θ ) listed in Table 1, one obtains θ (τp ) = I1 + I2 + αus (τp ) + β   = G us (τp ), I1 + I2 + αus (τp ) + β .

(29)

The pair of nonlinear equations governing the surface density and the subsurface layer concentration, Eq. (25) and a dimensionless dynamic isotherm listed in Table 1, are solved using a predictor–corrector method. The two integrals in Eq. (25) are evaluated numerically by partitioning the time domain τ into p parts of small intervals (τi = τi − τi−1 , i = 1, 2, . . . , p, and integrated separately by the trapezoidal rule over each interval [12]:

There is only one unknown, us (τp ), in the above equation and the root is found numerically by Newton’s method. Then, θ (τp ) is calculated by inserting us (τp ) into Eq. (28). The error incurred by numerical integration in determining the unknowns us (τp ) and θ (τp ) by the trapezoidal rule can be evaluated by expanding the last integral in Eq. (25) in a Taylor series from τp−1 to τp . Hence, the time truncation (p) error ET for the pth step is

p  2  F (τi ) + F (τi−1 ) √ √ θ (τp ) ≈ √ τi − τi−1 2 π

(p) ET ≡ θ¯ (τp ) − θ (τp )

i=1 p  us (τi )F (τp

− τi ) + us (τi−1 )F (τp − τi−1 ) 2 +√ 2 π i=1



× τp − τi − τp − τi−1 . (27) The summations in Eq. [27] can be further simplified as θ (τp ) = I1 + I2 + αus (τp ) + β,

(28)

2 ≈√ π

ζ =τp

us (ζ )F (τ − ζ ) d τp − ζ

ζ =τp−1

 2 us (τp ) + us (τp−1 )F ((τp ) − (τp −√ 2 π 3/2 

−(τp d  us (ζ )F (τ − ζ )
= √ (30) ζ =η 3 π dζ

314

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

where θ¯ (τp ) is the exact solution for the dimensionless surface density at time τp and τp−1
η
τp . (p) The corrector error ET incurred during time integration can be calculated by means of the predictor. In this paper, a (p) forward difference predictor is used to evaluate ET . The trapezoidal rule approximation is applied from 0 to τp−1 as before, but the last integral of Eq. (25) from τp−1 to τp is approximated by a forward difference. The resulting equation is

2 θF (τp ) = I1 + I2 − √ us (τp−1 )F ((τp ) (τp π = I1 + I2 + 2β,

(31)

where θF (τp ) is the predicted value of the surface density. (p) The error EF of the forward difference predictor is similarly obtained by expanding in a Taylor series the last integral of Eq. (25) from τp−1 to τp : (p) EF = θ¯ (τp ) − θF (τp )

2 =√ π

ζ =τp

us (ζ )F (τ − ζ ) d τp − ζ

ζ =τp−1



2 − √ us (τp−1 )F ((τp ) − (τp π 

2 d  =√ us (ζ )F (τ − ζ )
ζ =η π dζ

ζ =τp 3/2



d  us (ζ )F (τ − ζ )
ζ =η dζ

(p) = 4ET .

(p)

=

3/2 (τp+1 . 3/2 (τp

(τp

(32) (p)

qj

=

(34)

Hence, the value of the time-step size (τp+1 at the next step p+1 at the next can be chosen so that the corrector error ET

(p)

= f qj

(p−1)

+ (1 − f )qj

,

(37)

 ,

(38)

2 Nl2 ((ξj +1 + (ξj ) (p)  (p) uj +1 − uj ×

(ξj +1

−

(p)

(p)

uj − uj −1 (ξj

where (τp = τp − τp−1 , ξj is the location of the j th mesh (p) point, (ξj = ξj − ξj −1 , uj = u(τp , ξj ), j = 1, 2, . . . , Nx − 1, and ξ1 = 0 and ξNx = 1. The value of f in Eq. (37) depends on the time integration scheme: f = 1.0 for the backward difference scheme and f = 0.5 for the trapezoidal rule. The FDE (37) is second-order accurate with respect to spatial discretization. The initial conditions are θ (0) = 0,

(p)

ET

(p−1)

uj − uj

Comparing Eqs. (30) and (32) makes it plain that it is not necessary to know θ¯ (τp ), which cancels when one of these equations is subtracted from the other, to estimate the error, which is  1 (p) ET = θ (τp ) − θF (τp ) . (33) 3 The above algorithm of error estimation can also be used in choosing the time-step size adaptively [18]. From (p) 3/2 (p) 3/2 Eq. (33), ET is proportional to (τp , or ET = K(τp , where K is a parameter. Therefore, if K varies slightly over a very small time interval from τp to τp+1 such that it can be (p+1) expanded in a Taylor series in (τp , the value of ET can be estimated as (p+1) ET

The required accuracy in computed results to be reported in Section 3 is 10−3 or better. Since for the calculation of θ (τp ), one needs to determine the integrals in Eq. (25) from τ = 0 to τp , computational errors incurred during the first p − 1 steps would accumulate and affect the accuracy of the solution at the pth step. Thus, the suitable value of tolerance / for each step should be less than the required accuracy 10−3 , and needs to be investigated numerically.

After the space domain is discretized with Nx points, the finite difference approximation of Eq. (19) is generated by the control volume method [19]. The resulting finite difference equation (FDE) at the pth time step for the j th mesh point is

(ζ − τp−1 ) d τp − ζ

−4(τp √ = 3 π

In the calculations, a small time interval (τ1 = 10−4 is taken as the first step size. Subsequent step sizes are then chosen adaptively using Eq. (35). However, if a step size is too large or its value changes too much from one step to the next, the estimate provided by Eq. (34) can be inaccurate. Therefore, to keep the step size prediction more accurate, the selection of adaptive step sizes is further restricted by the following rules: (τp+1
2.0, (τmax
0.1. (36) (τp

2.3. Finite difference (FD) method

ζ =τ p−1

×

step is less than or equal to a specific tolerance /, or, from Eq. (34)   / 2/3 . (τp+1 = (τp (35) (p) ET

(0)

uj = 1,

j = 1, 2, . . . , Nx .

(39)

Because the above initial conditions entail a rapid change near ξ = 0, the first few time steps are integrated via the backward difference method to provide the necessary smoothing [20].

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

Since in the control volume method a midpoint approximation is used at each mesh point, a fictitious mesh point at ξ0 = −ξ2 is used to formulate the FDE for the interfacial boundary condition Eq. (21) at ξ = 0 that is also secondorder accurate with respect to spatial discretization, (p)

Nl

(p)

u − u0 θ (p) − θ (p−1) =f 2 (τp 2(ξ2

(p−1)

+ (1 − f )

u2

(p−1)

− u0 2(ξ2

,

(40)

where θ (p) ≡ θ (τp ). The FDE for the Dirichlet boundary condition Eq. (20) for point Nx at ξ = 1, at time τp , is (p)

uNx = 1.

(41)

The diffusion-controlled adsorption condition gives the remaining equation: 

(p) θ (p) = G u1 , θ (p) .

(42)

The resulting set of Nx + 2 nonlinear FDEs, Eqs. (37)–(42), are solved by Newton’s method. The choice of the location of the mesh points is important to both the accuracy and the efficiency of the calculations. Because the concentration has a large slope near ξ = 0 for small times τ , more mesh points are required near ξ = 0 to resolve the rapidly changing concentration profiles there. To resolve the sharp gradient, a “stretched” mesh is generated from a logically uniform grid in this paper as [21] ςi =

i −1 , Nx − 1

ξi =

eλςi − 1 , eλ − 1

i = 1, 2, . . . , Nx ,

(43)

where the parameter λ is a positive number. On physical grounds, it is appropriate to determine λ from the penetra√ tion depth ξd ≡ τ /Nl . At the first time step τ1 , the concentration has the largest gradient near ξ = 0, and is nearly a straight line over 0
ξ
ξd (τ1 ). Thus, the first few mesh points should be placed between ξ = 0 and ξ = ξd (τ1 ) so that the computed concentration profile is reasonably smooth near the interface. If n points are assigned to be within 0
ξ
ξd (τ1 ), then √ τ1 eλ(n−1)/(Nx −1) − 1 = . eλ − 1 Nl

error [18],

315

 + (τp 1 +

 (τp (p−1) q 2(τp−1 j  (τp (p−2) , − q (45) 2(τp−1 j   (p−1) (p−1) − u1 u2 (τp (τp (p) (p−1) θ¯ = θ + 1+ Nl 2(τp−1 (ξ2  (p−2) (p−2) − u1 (τp u2 − , (46) 2(τp−1 (ξ2

(p) u¯ j

(p−1) = uj

(p) where u¯ j and θ¯ (p) are the predicted values of the concentration at j th mesh point and surface density at time τp . The time truncation errors for the Adams–Bashforth predictor and the trapezoidal rule corrector at the pth step are evaluated by a Taylor series expansion analysis (detailed on p. 37 of Ref. [18]). Similarly to the procedure described in Eqs. (30)–(33), the time truncation errors for the corrector at the pth step are

θ (p) − θ¯ (p) , dp,0 ≡

(τp  3 1 + (τp−1

(p)

dp,j

(p)

uj − u¯ j . ≡

(τp  3 1 + (τp−1

(47)

With the error estimates for the pth time step above, the (p +1)th step size is chosen based on the fact that the relative error is proportional to (τ 3 . A maximum error of 0.1% per time step, / = 10−3 , is prescribed in the computations [18]. Therefore, the adaptive step size for the (p + 1)th time step is  1/3 / (τp+1 = (τp (48) , dp ∞ where dp ∞ ≡ maxj =0,...,Nx (|dp,j |). However, if the step size is too large or the ratio between two successive time steps is too big, the predicted step size can be inaccurate. Similarly to Section 2.2, to keep the algorithm for time-step size prediction accurate, the ratio between two successive time steps is restricted to be less or equal than 2.0, and the maximum step size is set to be (τmax = 0.1. 2.4. Finite element (FE) method

(44)

The value of n is chosen to be 3 in the calculations reported in Section 3, and λ in Eq. (44) is determined by Newton’s method. The error estimation at each step and the adaptive timestep size selection are based on the predictor–corrector method, as in the previous section. At first, four backward difference time steps with a fixed (τ = 10−4 are used to provide the necessary smoothing [20]. For subsequent time steps, a second-order Adams–Bashforth predictor is used with the trapezoidal rule corrector to evaluate the relative

The set of governing one-dimensional, time-dependent equations, Eqs. (19)–(24), is solved numerically with the method of lines using the Galerkin finite element method with quadratic basis functions for spatial discretization [22] and an adaptive finite difference method for time integration [23,24]. The spatial domain 0
ξ
1 is divided into NE elements. The concentration profile is then expanded in terms of a series of quadratic basis functions φ j (ξ ), u(ξ, τ ) =

Nx  j =1

uj (τ )φ j (ξ ),

(49)

316

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

where uj s are unknown coefficients to be determined and Nx = 2NE + 1 is the total number of nodes. The Galerkin weighted residuals of Eq. (19) are constructed by weighting Eq. (19) with the quadratic basis functions φ i (ξ ), integrating the resulting expressions over the computational domain, and rearranging the integrals by integration by parts:  1  ∂u ∂ 2 u i − 2 φ dξ R = Nl2 ∂τ ∂ξ i

0

 1  ∂u dφ i 2 ∂u i φ + dξ Nl = ∂τ ∂ξ dξ 0

− φi


∂u

ξ =1 = 0, ∂ξ
ξ =0

i = 1, 2, . . . , Nx .

3. Results and discussion (50)

The weighted residual equation, Eq. (50), can be simplified by means of the boundary conditions Eq. (20) and (21): Ri =

 1  ∂u ∂u dφ i Nl2 φ i + dξ ∂τ ∂ξ dξ 0

dθ = 0, dτ − 1 = 0.

− δi1 Nl R Nx = uNx

i = 1, . . . , Nx − 1, (51)

An additional residual equation arises from a dimensionless isotherm listed in Table 1, 

(p) R 0 = θ (p) − G u1 , θ (p) = 0. (52) The initial conditions are θ (0) = 0,

u(0) j = 1,

j = 1, 2, . . . , Nx .

(53)

The FE grid is generated in a manner that parallels that in which the FD grid is generated. First, the size of each element is determined by the stretched grid generation method and then the middle node of each element is located halfway between the two end nodes, viz., ςi =

i , NE

ξ1 = 0,

i = 1, 2, . . . , NE , ξ2i+1 =

used to provide the necessary smoothing [20]. In the subsequent steps, a second-order Adams–Bashforth predictor is used with the trapezoidal rule to evaluate the relative time truncation errors and time-step sizes adaptively. The algorithm for choosing time-step sizes adaptively and step-size restrictions is identical to that detailed in Section 2.3. Finally, with the FE method with the stretched grid and the adaptive time integration in place, the system of Nx + 1 nonlinear algebraic equations in Eqs. (51) and (52) is solved by Newton’s method.

eλςi − 1 , eλ − 1

ξ2i−1 + ξ2i+1 . (54) 2 The parameter λ is determined by placing one element within the region 0
ξ
ξd (τ1 ): √ τ1 eλ/NE − 1 = . (55) λ e −1 Nl

ξ2i =

The time derivative of the concentration u at node j and that of the surface density are formulated with either the backward difference method or the trapezoidal rule over a time interval (τp = τp − τp−1 . Initially, four backward difference time steps with fixed time-step size (τ = 10−4 are

3.1. Code validation with the Henry isotherm The codes implementing the algorithms based on the I, FD, and FE methods have been written in FORTRAN. The calculations have been carried out on a Sun Ultra 60 workstation at Purdue University. For code validation and studying the effects of time-step sizes and the number of mesh points on solution accuracy, the codes have been tested for the limiting case of a linear, or Henry, isotherm for which an exact solution exists [16], ∞ 

e−βn τ/Nl θ (τ ) = 1 − 2Nl , β 2 + Nl + Nl2 n=1 n 2

2

(56)

where the βn ’s are the roots of the equation βn tan(βn ) = Nl . The convergence of different methods is tested by comparing the predicted values of the dimensionless time τx for dimensionless surface density θ to be x% saturated with that obtained from Eq. (56). Three points, 25%, 50%, and 95%, are of particular interest. A solution is said to have converged if the third decimal digit of τx stops changing upon further reduction in time-step size or increase in number of mesh points used. Table 2 summarizes results obtained for the case of the dynamic Henry isotherm with Nl = 1.0. The errors in the I method solutions are due to time truncation errors that arise from time integration. Thus, the accuracy of the I method solutions with fixed step sizes is improved by decreasing (τ . A step size (τ of 0.01 or less is required to obtain a precision of three digits after the decimal point. Similarly, for the adaptive I method solutions, decreasing the prescribed tolerance /, which represents the time truncation error for each step, increases the accuracy; / should be 10−5 or smaller to achieve three-decimal-digit precision. For numerical solutions obtained with the FD and FE methods, time truncation errors are kept below a prescribed tolerance by adaptive time integration. Hence, overall errors arise mainly from spatial discretization error and solution accuracy can be improved by using more grid points in the calculations. Table 2 shows that as the number of mesh points in the FD method is increased, the solutions become more accurate and finally converge. Table 2 makes it plain

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

Table 2 Convergence tests for the Henry isotherm with Nl = 1.0: time τx for θ to be x% saturated

τ25

τ50

τ95

0.077

0.512

3.622

IM with fixed step size (τ 0.050 0.020 0.010 0.005

τ25 0.073 0.076 0.077 0.077

τ50 0.508 0.511 0.512 0.512

Table 3 Computational requirements and code execution times with the Henry isotherma Method (Nl )

Exact solutions (Eq. (56))

IM with adaptive step size τ95

/

τ25

τ50

τ95

3.624 3.622 3.622 3.622

10−3

10−4 10−5

0.077 0.077 0.077

0.510 0.512 0.512

3.633 3.626 3.622

FDM with uniform grid

FDM with stretched grid

Nx

τ25

τ50

τ95

Nx

τ25

τ50

τ95

40 80 160 320

0.077 0.077 0.077 0.077

0.513 0.513 0.512 0.512

3.621 3.621 3.622 3.622

10 20 40 80

0.074 0.077 0.077 0.077

0.497 0.510 0.512 0.512

3.655 3.626 3.622 3.622

FEM with uniform elements

FEM with stretched elements

NE a

τ25

τ50

τ95

NE a

τ25

τ50

τ95

10 20 40 80

0.078 0.077 0.077 0.077

0.513 0.513 0.512 0.512

3.621 3.621 3.621 3.622

2 5 10 20

0.079 0.077 0.077 0.077

0.514 0.512 0.512 0.512

3.621 3.622 3.622 3.622

a N = 2N + 1. x E

that, in this example, 160 mesh points when a uniform mesh is used, and only 40 mesh points when a stretched mesh is used, suffice for FD solutions so that all three values of τx stop changing beyond the third digit after the decimal point. Table 2 shows that similar results are obtained with the FE method; however, convergence is achieved with fewer mesh points with the FE method than with the FD method. Table 2 makes it plain that in this example, 40 elements, or equivalently 81 mesh points, when a uniform mesh is used, and 5 elements, or 11 mesh points, when a stretched mesh is used, suffice with the FE method. Table 3 reports similar convergence tests that have been done to evaluate the performance of the three methods for different values of the diffusion layer thickness Nl . Since the required accuracy is 10−3 , the calculations are stopped when θ = 0.999. Table 3 also reports CPU times texe in central processor time units required by the calculations. When analyzing dynamic adsorption with the I method, the adaptive version of the method can guide wisely finding of appropriate step sizes to save computation time while ensuring accuracy. Table 3 shows that when Nl = 0.1 and 1.0, because adsorption is fast and hence only some hundreds of time steps are needed to reach equilibrium, solutions obtained using the I method with either fixed or adaptive time-step sizes require roughly the same amount of computer time. However, when Nl = 10, because adsorption is slow and hence many more time steps are needed by the version of the code using fixed step sizes, about 12 min are required when the fixed step size version of the code is used, whereas only 15 s are needed by the version that determines step sizes

317

IM with fixed step size (τmax texe , s IM with adaptive step size /max texe , s FDM with uniform grid Nx, min texe , s FDM with stretched grid Nx, min texe , s FEM with uniform grid NE, min b texe , s FEM with stretched grid NE, min b texe , s

0.1

1.0

10

10−3 2.64

10−2 1.55

10−2 748.60

10−5 3.20

10−5 1.75

10−5 15.04

40 0.02

160 0.20

2560 59.13

10 0.02

40 0.08

160 3.30

10 0.02

40 0.17

160 10.45

5 0.01

5 0.04

20 0.69

a t : Code execution time in seconds for computation to reach exe equilibrium (θ = 0.999). b N x, min = 2NE, min + 1: Grid points needed for the convergence at the third digit after decimal point.

adaptively. Hence, since the adaptive I method algorithm is more efficient, it is adopted in the calculations that use the Langmuir and Frumkin isotherms to be reported in the next section. Table 3 makes it plain that for the FD and FE methods, the algorithms using stretched grids are more efficient than ones using uniform grids, because the former require fewer grid points and consequently consume less computer time. The benefits of using stretched grids are more recognizable, especially for the larger value of Nl . For example, in the situation where Nl = 10, whereas the FD algorithm using a stretched grid requires only about 3 s, that using a uniform grid requires about 20 times longer or 60 s. Therefore, for reasons of computational efficiency, results to be reported in the next section, where the isotherms are nonlinear, have been attained with FD and FE algorithms using stretched grids. Among the three methods, solutions obtained with the I method even when time-step sizes are determined adaptively have always been found to require more computer time than solutions obtained with the FD or FE methods using stretched grids. Furthermore, algorithms based on the FD or FE methods determine concentration profiles simultaneously with surface densities. Therefore, codes based on the FD and FE methods are not only more efficient but more informative than codes based on the I methods. 3.2. Solutions for the Langmuir isotherm The convergence and execution times required by the three methods are evaluated in this section when the Lang-

318

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

Table 4 Computational requirements and code execution times with the Langmuir isotherm Nl

NcL

τ95

0.1 0.1 0.1 1 1 1 10 10 10

1 10 100 1 10 100 1 10 100

0.202 0.114 0.094 2.401 1.108 0.675 31.318 2.709 0.819

IM

FDM

FEM

/max

texe , s

Nx, min

texe , s

NE, min

Nx, min

texe , s

10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5

1.21 0.32 0.06 1.48 0.96 0.43 10.76 4.95 1.50

10 10 10 40 40 40 160 160 160

0.04 0.04 0.04 0.08 0.07 0.07 3.45 2.33 0.88

5 5 5 5 5 5 20 20 20

11 11 11 11 11 11 41 41 41

0.01 0.01 0.01 0.03 0.02 0.03 1.05 0.63 0.76

Adaptive time integration and stretched grids are used.

muir isotherm is used. Table 4 lists the requirements of the algorithms based on the I, FD, and FE methods to attain convergence in τx of three digits after the decimal point. Since there is no analytical solution to the problem when the Langmuir isotherm is used, the accuracy of the solutions is tested not only by convergence tests but also by requiring that the values of τ25 , τ50 , and τ95 calculated by the three methods all converge to the same numbers. For brevity, Table 4 only lists values of τ95 . Since for NcL  1, the solutions with the Langmuir isotherm are very close to those with Henry’s isotherm, the convergence of solutions with the Langmuir isotherm has been tested only over the range of 1
NcL
100. Moreover, the range of diffusion layer thicknesses considered vary as 0.1
Nl
10. For solutions obtained from the I method, the prescribed error /max has always been found to be 10−5 for achieving 0.1% accuracy, as shown in Table 4. For solutions obtained with the FD and FE methods, the required minimum number of grid points has been found to be the same for all NcL values, but to increase with increasing values of Nl . Moreover, for the same value of Nl , the required value of /max or the minimum number of grid points for the Langmuir isotherm has been found to be about the same as that for Henry’s isotherm, implying that the calculation requirements are affected primarily by Nl . Figure 1a shows the calculated dimensionless surface densities determined by the I, FD, and FE methods. The larger the value of NcL , the faster is the adsorption. Solutions obtained using the I, FD, and FE methods are so close that they cannot be distinguished visually from one another in Fig. 1a. However, the code execution times for the three methods are dramatically different. To attain the same accuracy, the FE method is the most efficient method, as shown in Table 4. When adsorption is fast (larger NcL or smaller Nl ), the difference of different methods in execution times is quite small (less than 1 s). But when adsorption is slow, the code based on the FE method is much more efficient than those based on the I and FD methods. For example, in the case of Nl = 10 and NcL = 1, about 1 s is needed by the FE algorithm, 3 s by the FD algorithm, and 11 s by the algorithm based on the I method.

Fig. 1. Evolution in time of calculated (a) surface densities and (b) surface pressures when the Langmuir isotherm is used for the case in which Nl = 1: (1) NcL = 100; (2) NcL = 10; (3) NcL = 1 and (4) NcL = 0.1. Solution methods are: adaptive I method with / = 10−5 ; FD method with a stretched grid and Nx = 40; FE method with a stretched grid and NE = 5 or Nx = 11. The small differences between solutions obtained with the three methods are not visually recognizable from the figures.

In practice, it is often the case that surface pressures are measured and of more interest than surface densities. Diffusion coefficients are either obtained from independent measurements or chosen by the best fit of surface tension data to computed results [10]. The dimensionless surface pressure Π can be calculated directly from the surface density θ , as shown in Table 1. The times to attain 50% and 95% reduction in surface pressure, i.e., τΠ50 and τΠ95 , are

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

Fig. 2. Effect of Nl on the evolution in time of surface density when the Langmuir isotherm is used. Here NcL = 10 and θ is calculated by the FE method for different values of the diffusion layer thickness: (1) Nl = 0.1; (2) Nl = 1; (3) Nl = 10; and (4) Nl = 100.

the most commonly used indices. Thus, the convergence of τΠ50 and τΠ95 has also been studied. The requirements for accuracy in three decimal digits in the I, FD, and FE methods are the same as those for θ listed in Table 4. Figure 1b shows surface pressures calculated from the various methods. Once again, the three methods give visually the same solutions and the time to reach equilibrium surface tension decreases as the value of NcL increases. Figure 2 illustrates the effect of Nl on the evolution in time of surface density. When Nl
10, the rate of adsorption increases as Nl decreases. The equilibration time (τ99 ) when Nl = 1.0 (curve 2) is about 10 times larger than that when Nl = 0.1 (curve 1). However, when Nl  10, an increase in Nl has a negligible effect on the solutions. In fact, the solutions when Nl = 10 and Nl = 100 overlap, as can be seen from curves 3 and 4 in Fig. 2. Therefore, in this example, the interfacial adsorption is the same for l = 10L and l = 100L. Figure 3 further explores the effect of Nl on adsorption by comparing the equilibration times for several values of the bulk concentration NcL . Since surface densities are not easily observed directly in experiments, surface tension is a more reliable measure for detecting the adsorption. The equilibration time for adsorption at a clean surface is taken as that at which a 95% reduction in surface tension has occurred. Thus, τΠ95 is adopted as the equilibration time τe in this paper. Figure 3 shows the effect of NcL on τe . In the case of Nl  1, e.g., Nl = 1000, increasing the value of NcL accelerates the adsorption. This occurs, as pointed out by Lin et al. [25], because at the same dimensionless surface density θ , the larger NcL is, the smaller is the sublayer concentration us given by the dimensionless Langmuir isotherm (Table 1). Hence, the larger NcL is, the larger is the diffusive flux for any value of θ and consequently the smaller is the equilibration time. For the case of Nl = 0.1, which represents a small diffusion layer thickness, τe starts from 0.3 in the Henry region (NcL
0.1), and decreases slightly and tends to 0.1 as NcL increases. With increasing Nl , the adsorption process takes longer to equilibrate. It is not until

319

Fig. 3. Variation of equilibration time τe for the Langmuir isotherm with concentration NcL : solid line, Nl = 100 and 1000 (cannot be distinguished); dashed line, Nl = 10; dash–dot line, Nl = 1; and dotted line, Nl = 0.1. Table 5 Computational requirements and code execution times with the Frumkin isotherm when Nl = 10 and NcF = 10 A = −3

A=0

A=3

IM (τmax texe , s τ95

10−5 0.73 0.786

10−5 4.95 2.709

10−5 7.63 8.436

FDM Nx, min texe , s τ95

160 0.50 0.786

160 2.33 2.709

160 2.99 8.436

FEM NE, min Nx, min texe , s τ95

20 41 0.42 0.786

20 41 0.63 2.709

20 41 0.86 8.436

Adaptive time integration and stretched grids are used.

Nl = 100 that the τe vs NcL curves in the low-concentration region (NcL
1.0) stop changing with increasing Nl . In the higher concentration region (NcL  10), the equilibration times τe for both Nl = 10 and Nl = 100 coincide. Therefore, for dynamic adsorption with NcL  10, Nl = 10 is a good enough approximation of infinite diffusion length. 3.3. Solutions for the Frumkin isotherm The Frumkin isotherm accounts for certain binary molecular interactions between adsorbed surfactant molecules. For many surfactants, it is capable of describing the time variation of dynamic surface tension more accurately than the Langmuir isotherm [26]. When the binary interaction parameter A = 0, in other words when there are no net energetic interactions between surface molecules, the Frumkin isotherm reduces to the Langmuir isotherm. Table 5 illustrates the requirements of the three numerical methods for attaining convergence in τx when the Frumkin isotherm is used. At the same values of Nl , the requirements

320

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

calculated surface tensions show larger variations, as shown in Fig. 4b. The onset of tension reduction takes off more sharply in the attractive case, A = −3, than in the other two cases. Furthermore, the time for the surface tension to reach 95% reduction is almost the same in all three cases. Therefore, the effect of molecular interactions is more pronounced at the beginning of surface tension reduction rather than for the times approaching the equilibrium. 3.4. Possible extensions to other models The aforementioned methods can also be applied, with proper modifications, to the following problems and help improve previously published solution methods. (a) The mixed kinetics problem in planar coordinates in one dimension [11] can be readily analyzed by incorporating simple and straightforward changes into the three algorithms presented here. (b) The algorithms presented in this paper can readily be extended to analyze both the diffusion-controlled or the mixed kinetic problem in spherical coordinates in one dimension [13], as in problems involving drops or bubbles. (c) The FD and FE algorithms to analyze the problems described in (b) can also be readily adapted to account for pulsating area conditions [15]. Fig. 4. Evolution in time of calculated (a) surface densities and (b) surface pressures when the Frumkin isotherm is used for the case in which Nl = 10 and NcF = 10: solid line, A = 0; dash–dot line, A = −3; and dashed line, A = 3.

for different values of A are about the same as those for the Langmuir isotherm (A = 0). Since there is no analytical solution to the problem when the Frumkin isotherm is used, accuracy is tested by requiring that τ25 , τ50 , and τ95 from various methods converge to the same answers. Table 5 reports only values of τ95 . The code execution times for various methods given in Table 5 show that the FE code is the most efficient. The FE code takes less than 1 second to solve the problem with the Frumkin isotherm for the parameters given in Table 5. Figure 4 demonstrates the effect of the molecular interaction parameter A on the evolution in time of θ and Π . Once again, the codes based on the I, FD, and FE methods give almost the same solutions when the Frumkin isotherm is used. Therefore, since the solutions obtained from the various methods are indistinguishable from one another, only the FE solutions are shown in Fig. 4. When A > 0, i.e., when the surfactant molecules exhibit net repulsive interactions, adsorption is slower than when A = 0. In contrast, when A < 0, i.e., when the surfactants molecules exhibit net attractive interactions, the adsorption is faster than when A = 0. Initially (θ
0.5), the calculated surface densities are almost the same, as shown in Fig. 4a. When θ  0.5, the influence of the A parameter can be seen. However, the

For more complicated two-dimensional problems involving pulsating drops or bubbles, we feel that the FE method is expected to provide the most accurate solutions based on recent experience with problems involving drop breakup [27,28].

4. Conclusions Three methods for solving the problem of one-dimensional diffusion-controlled dynamic adsorption in planar coordinates are analyzed. A technique for error estimation in the integral (I) method with the trapezoidal rule approximation is developed. Error estimation is then used to choose adaptively time-step sizes in the I method. The analytical solution of the linear Henry isotherm, used for code validation and convergence tests, shows that the I method with adaptive time-step sizes is more efficient to achieve the same accuracy, with respect to code execution time, than that with uniform time-step sizes. The maximum tolerance / for the I method at each step to achieve 0.1% accuracy is determined to be 10−5 . The FD and FE algorithms that are developed use a predictor–corrector method with adaptive time stepping to control time truncation error and a stretched grid generation technique to minimize spatial discretization error. Convergence tests with Henry’s isotherm show that much fewer grid points are needed to achieve an accuracy of 0.1% using a stretched grid than a uniform grid. Once the I, FD, and FE

Y.-C. Liao et al. / Journal of Colloid and Interface Science 258 (2003) 310–321

methods are developed, analyzed, and optimized, they are used in the solution of the problems with the Langmuir and Frumkin isotherms, which are nonlinear, as examples. All three methods converge to the same solution, which cannot be obtained analytically. The effects of the dimensionless finite diffusion layer thickness Nl and the dimensionless concentration Nc are also examined. Efficiency of these three methods is also compared with respect to code execution times. At high concentrations or under fast adsorption conditions, all three methods require almost the same amount of computer time. However, when adsorption is slow, the FE method is the most efficient.

Acknowledgments This research was supported in part by grants from the National Science Foundation (Grant CTS 96-15649, CTS 96-15649, and 0135317) and the National Institutes of Health (Grant HL 54641-02) to Elias I. Franses, and a grant from the Basic Energy Sciences Program of the US Department of Energy (Grant DE-FG02-96ER14641) to Osman A. Basaran.

References [1] L. Brown, G. Narsimhan, P.C. Wankat, Biotechnol. Bioeng. 36 (1990) 947. [2] J.E. Valentini, W.R. Thomas, P. Sevenhuysen, T.S. Jiang, H.O. Lee, L. Yi, S.-C. Yen, Ind. Eng. Chem. Res. 30 (1991) 453.

321

[3] R.H. Notter, P.E. Morrow, Ann. Biomed. Eng. 3 (1975) 119. [4] B. Ambravaneswaran, O.A. Basaran, Phys. Fluids 11 (1999) 997. [5] R. Defay, G. Petré, in: E. Matijevi´c (Ed.), Surface and Colloid Science, Vol. 3, Wiley, New York, 1971. [6] A.F.H. Ward, L. Tordai, J. Chem. Phys. 14 (1946) 453. [7] K.L. Sutherland, Aust. J. Sci. Res. Ser. A 5 (1952) 683. [8] R. Miller, Colloid Polym. Sci. 259 (1981) 375. [9] B.J. McCoy, Colloid Polym. Sci. 261 (1983) 535. [10] C.-H. Chang, E.I. Franses, Colloids Surf. A 100 (1995) 1. [11] C.-H. Chang, E.I. Franses, Colloids Surf. 62 (1992) 321. [12] R. Miller, G. Kretzschmar, Colloid Polym. Sci. 258 (1980) 85. [13] S.-Y. Lin, K. McKeigue, C. Maldarelli, AIChE J. 36 (12) (1990) 1785. [14] K. Wantke, K. Malysa, K. Lunkenheimer, Colloids Surf. A 82 (1994) 183. [15] C.-H. Chang, E.I. Franses, Chem. Eng. Sci. 49 (3) (1994) 313. [16] K.J. Mysels, H.L. Frisch, J. Colloid Interface Sci. 99 (1984) 136. [17] J.K. Ferri, K.J. Stebe, Adv. Colloid Interface Sci. 85 (2000) 61. [18] P.M. Gresho, R.L. Lee, R.C. Sani, in: C. Taylor, K. Morgan (Eds.), Recent Advances in Numerical Methods in Fluids, Vol. 1, Pineridge, Swansea, 1979, pp. 27–79. [19] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw– Hill, New York, 1980. [20] M. Luskin, R. Rannacher, Applicable Anal. 14 (1982) 117. [21] P. Knupp, S. Steinberg, The Fundamentals of Grid Generation, CRC Press, Boca Raton, FL, 1993. [22] G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice–Hall, Englewood Cliffs, NJ, 1973. [23] O.A. Basaran, P.M. Burban, S.R. Auvil, Ind. Eng. Chem. Res. 28 (1) (1989) 108. [24] X. Zhang, R.S. Padgett, O.A. Basaran, J. Fluid Mech. 329 (1996) 207. [25] S.-Y. Lin, K. McKeigue, C. Maldarelli, Langmuir 10 (1994) 3442. [26] C.-T. Hsu, C.-H. Chang, S.-Y. Lin, Langmuir 16 (2000) 1211. [27] E.D. Wilkes, S.D. Phillips, O.A. Basaran, Phys. Fluids 11 (1999) 3577. [28] P.K. Notz, A.U. Chen, O.A. Basaran, Phys. Fluids 13 (2001) 549.