Mass transfer coefficients for laminar longitudinal flow in hollow-fibre contactors

We consider the problem of predicting the rate of mass transfer to a fluid flowing parallel to the axes of randomly placed aligned tubes, a model of hollow-fibre contactors. The analysis is carried out for the limiting cases of short contactors, for which the concentration boundary layers remain thin compared with the radius of the tubes, and for the fully developed case corresponding to very long tubes. Numerical simulations for random arrays are carried out for $N$ randomly placed tubes within a unit cell of a periodic array. It is shown that the mass transfer coefficient for the fully developed case is vanishingly small in the limit $N\rightarrow \infty$. This suggests that the mass transfer coefficient for a random array of tubes of radius $a$ enclosed in a shell of radius $S$ will vanish logarithmically as the ratio $S/a$ is increased. This behaviour arises due to the logarithmically divergent nature of concentration disturbances caused by each tube in the plane normal to its axis. A theory is developed for determining conditionally averaged velocity and concentration fields and its predictions are shown to compare very well with the results of rigorous numerical computations. The predictions of the theory are also shown to compare well with the measurements of the mass transfer coefficients in hollow-fibre contactors reported in the literature.


Introduction
Rapid progress in membrane technology over the last few decades has made it possible to replace many of the conventional separation processes, such as gas absorption, distillation, and liquid-liquid extraction, by more efficient membranebased processes. Typical equipment consists of a module of hollow fibres made out of a suitable membrane material placed inside a cylindrical vessel or shell with a gas or liquid flowing on either side of the membrane as shown in figure 1. The difference in the concentration of solute on the two sides causes solute to diffuse from the bulk of the fluid on one side to that on the other through the membrane, thus effecting the separation. The entrance and the exit of the shell-side fluid may be located such that the bulk of the fluid flow on the shell side is either perpendicular to the tube bundle (cross-flow) or parallel (longitudinal). Such hollow-fibre contactors are now commonly used, for example, in dialysis, gas separation, and blood oxygenators. The mass transfer process is often controlled by the resistance on the shell side and accurate estimates of the shell-side mass transfer coefficients are needed for designing these units.
Although shell-and-tube configurations have been widely used in heat exchangers for over a century, the correlations typically used for designing heat exchangers are not very useful in designing hollow-fibre contactors. This is so in spite of the similarity in the equations governing the heat and mass transfer processes because of a few important differences: (i) the hollow-fibre modules typically operate at Reynolds numbers that are small compared to those encountered in heat exchangers; (ii) the baffles placed on the shell side to support tubes and create significant mixing and cross-flow in heat exchangers are generally absent in the fibre contactors; and (iii) whereas the tubes are arranged in a regular, square or hexagonal, array in most heat exchangers, the fibre bundles are typically randomly arranged. The fibres are also flexible and often not parallel to each other. These differences have prompted several investigators to measure mass transfer coefficients in hollow-fibre modules in recent years (Yang & Cussler 1986;Prasad & Sirkar 1988;Costello et al. 1993). Theories for predicting mass transfer coefficients have also been developed, with the most notable contribution being made by Bao, Liu & Lipscomb (1999) who determined the mass transfer coefficients for short modules for which the concentration boundary layers are very thin compared to the tube radius. The tubes were assumed to be rigid and parallel to each other and the flow was assumed to be fully developed, laminar and parallel to the tubes. For this case the mass transfer coefficient can be determined from knowledge of the wall stress distribution, which they determined numerically for random arrays. Unfortunately, the predictions of this theory do not compare well with the values measured by the investigators cited above. One reason for this might be the limitation of the theory to short contactors since the experiments were carried out for systems of practical interest in which the length of the fibres is often much greater than the radius. It should be noted that the fibres may be arbitrarily close to each other in random arrays so that the assumption of non-overlapping concentration boundary layers might be too restrictive from the practical point of view. The overlap will reduce the mass transfer coefficient and one might therefore expect the theory to overpredict the mass transfer rates.
The correlations for the mass transfer coefficient obtained by the investigators cited earlier appear to suggest that the mass transfer coefficient is proportional to Re n with the exponent n being in the range 0.6-0.9 for the Reynolds number Re in the range 0-300 whereas the theory due to Bao et al. (1999) gives n = 1/3. Thus one might be tempted to discard the effect of overlapping boundary layers and instead suggest that more complicated flow models are necessary to explain the measured values of the mass transfer coefficients. Indeed it has been suggested in the literature that flow channelling, turbulence, flexibility of the fibres, and the presence of cross-flow near the entrance and exit of the contactors may be responsible for the more complicated observed behaviour of the mass transfer coefficients (see, for example, Yang & Cussler 1986;Costello et al. 1993;Bao et al. 1999;Wu & Chen 2000). A careful inspection of the data, however, reveals that the non-dimensional mass transfer coefficients measured by these investigators are too small for turbulence effects to be at work. Indeed, at least in one case, an extremely low mass transfer rate was measured by Yang & Cussler (1986). It is possible that at lower Reynolds numbers a very small fraction of the overall length of the contactor has boundary layers that are thin, resulting in a much smaller mass transfer coefficient than that at larger Reynolds numbers where a greater fraction of the fibre length contributes to the mass transfer. It is important therefore to carry out a systematic analysis for finite-length contactors based on a simple model of fully developed laminar flow before more complicated flow models are invoked. This is the purpose of the present study. It will be shown that the calculations based on laminar flow can yield estimates that are in good agreement with the measured values.
Rigorous calculations of the shell-side mass transfer coefficient in finite-length contactors are quite difficult even for the simple case of laminar flow parallel to a random array of straight rigid tubes. Therefore we shall analyse in detail first the two limiting cases of short contactors and long-contactors. For short contactors the analysis will be similar to that of Bao et al. (1999) but will account for the overlap of the concentration boundary layers in a somewhat ad hoc manner. The long-contactor calculations turn out to be more involved and interesting from the physics and theory standpoint so that the major portion of this paper is devoted to this limit.
For the long-contactor case we examine in detail a special situation in which the volumetric flow rates of the shell-and tube-side fluids are equal. For this case the average concentration profiles on the tube and shell sides are linear and the mass transfer across the tubes is independent of the axial position, a case equivalent to the constant wall heat flux examined by Sparrow, Loeffler & Hubbard (1961). We consider the case of random arrays by generating hard-disk configurations of N non-overlapping disks placed within a two-dimensional unit cell of a periodic array and computing mass transfer coefficient for each configuration. We find that the mass transfer coefficient vanishes logarithmically in the limit N → ∞, suggesting that the mass transfer coefficient for a truly random array is vanishingly small. This surprising logarithmic divergence arises from the fact that the fundamental singularity for the Laplace equation in a two-dimensional space corresponding to a point source is log r, r being the distance from the source. For the fully developed case the concentration inside and outside the tubes satisfies the Poisson equation with the strength of the source or sink related to the fluid velocity. The conditionally averaged concentration field obtained by averaging over all configurations of tubes with the position of one of the tubes fixed satisfies the Laplace equation at great distances from the tube. Detailed analysis of the source and sink near the tube, however, shows a non-zero net source leading to the logarithmic divergence. For a finite number of tubes in a periodic array this net source is balanced by a sink uniformly distributed over the unit cell, leading to a mass transfer coefficient that scales inversely with the logarithm of the unit cell size. For a random array bounded by a shell of radius S, the mass transfer coefficient will vanish logarithmically with the increase in the ratio S/a, a being the outer radius of the fibres.
The result that large fibre modules perform poorly does not seem to have been pointed out in the literature since there appears to be only one set of careful measurements of Sherwood number for the fully developed case with a large number of fibres. Yang & Cussler (1986) carried out experiments with 2100 fibres and found the Sherwood number, the non-dimensional mass transfer coefficient, to equal only 0.08. They were somewhat surprised by the result and hypothesized that channelling or uneven flow distribution on the shell side was responsible for such a low Sherwood number. We believe that this was not the case since our calculations give essentially the same estimate of the Sherwood number for their system. An exact expression is derived for the coefficient of logarithmic dependence that depends on the conditionally averaged velocity outside the fixed tube. An effectivemedium approximation is used for predicting the conditionally averaged velocity and concentration fields, and hence the coefficient of logarithmic dependence of the inverse of Sherwood number. The theory is shown to be in good agreement with the simulation results.
This paper is organized as follows. The governing equations and the method are described in § 2. The theory and results for random arrays are presented in § 3. The case of mass transfer at small axial distances is presented in § 4. The theory is compared with the experiments in § 5 and a summary of the work is given in § 6.

Formulation of the problem and the method
We begin with the analysis for the mass transfer coefficient at large axial distances. We consider a countercurrent shell-and-tube configuration as shown in figure 1. The flow on both sides is assumed to be laminar, fully developed, and parallel to the axes of the rigid tubes. We shall consider here only the special case when the volumetric flow rate of the shell-side fluid equals that on the tube side. The average concentrations of the solute on the tube and shell sides increase linearly with the axial distance in this case and the total mass transfer per unit length is independent of the axial position. With no loss of generality, we let the average concentration gradient equal 1/P e and write Here, the x 3 -axis is taken to be along the axes of the tubes and (x 1 , x 2 ) are the coordinates of a point in the plane normal to the tubes. The distances are nondimensionalized by a, the radius of the tubes. P e = aU/D s is the Péclet number based on the flow outside the tubes, where U is the superficial velocity of the fluid on the shell side, and D s is the diffusivity of the solute in the shell-side fluid. Upon substitution of (2.2) into the mass conservation equation for the solute we obtain for the shell side 3) where u s is the velocity of the fluid non-dimensionalized by the superficial velocity U and is the Laplacian operator in the (x 1 , x 2 )-plane. Since we have taken the volumetric flow rates of the shell-and tube-side fluids to be equal, the average velocity U t of the fluid inside a tube equals U/φ, where φ is the area fraction of the tubes. The mass conservation equation for the solute in the tube-side fluid reduces to where α c = D t /D s is the ratio of diffusivities in the two fluids, and u t is the standard non-dimensional parabolic profile for laminar flow through circular tubes. For a tube centred at the origin, we have Here, r is the radial distance from the centre of the tube. The negative sign on the right-hand side of (2.5) accounts for the countercurrent nature of the flows on the tube side. The positions of the centres of N tubes will be denoted by x α , α = 1, 2, . . . , N. These centres lie within a unit cell of a periodic array. The boundary conditions for the concentration of the solute are therefore spatial periodicity and the continuity of concentration and flux at the surface of the tubes: Note that we have assumed that the tube membrane thickness is negligibly small and offers no resistance to the mass transfer. We have also assumed that the partition coefficient of the solute in the two fluids is unity. The problem formulated here is similar to the heat transfer problem formulated by Sparrow et al. (1961). These investigators considered the special case of a periodic array (N = 1) and constant wall heat flux at the tube walls for which it is not necessary to solve for the temperature inside the tube. Equations (2.5) and (2.7) must be discarded and instead a boundary condition of uniform heat (or mass) flux at tube wall must be imposed in this case. The physics of the two problems are very similar and the heat or mass transfer coefficients for the two problems will have a similar dependence on the number of tubes.
The more general case of unequal volumetric flow rates is somewhat more complicated to analyse but it can be shown that the results obtained for the special case considered here will still be qualitatively applicable (Koo 2002).
We shall be interested in the Sherwood number, the non-dimensional mass transfer coefficient. The overall Sherwood number is defined as where Q is the rate of mass transfer per tube per unit length of the contactor, h ov is the mass transfer coefficient, and C ov is the difference between the average solute concentrations on the tube and shell sides. We shall use two kinds of averages. The first is a spatial average and the second is a fluid-velocity-weighted average, referred to in the literature as the mixing-cup average, Here, τ is the area of the unit cell non-dimensionalized by a 2 and A s is the area occupied by the shell-side fluid. Note that the integral in (2.10) is divided by τ only since the velocity u s is non-dimensionalized by the superficial velocity so that A s u s dA = τ. (2.11) The average concentrations for the tube-side fluid are defined in a similar manner.
The mass transfer per tube can be related to the average concentration gradient. Thus, it is easy to show that (2.12) Substituting for Q in (2.8) we obtain (2.13) 2.1. Determination of the velocity field We shall use the method of multipole expansion for determining the velocity and concentration fields. The method uses periodic fundamental singular solutions of the Laplace and biharmonic equations and their derivatives to construct velocity and concentration fields. We shall describe in more detail the procedure for determining the velocity field here, which follows the analysis presented in Sangani & Yao (1988b).
The shell-side fluid velocity satisfies where G is the pressure gradient non-dimensionalized by µU/a 2 . A multipole expansion expression for the velocity field is given by (Sangani & Yao 1988b) where A α n andÃ α n are the 2 n -multipoles induced by the presence of tube α,Ã 0 ≡ 0, and ∂ n k = (∂ n /∂x n k ) (k = 1, 2) is a short-hand notation for the nth-order partial derivative with respect to x k . The function S 1 is a spatially periodic function satisfying (Hasimoto 1959) (2.16) In the above expression, the x L are the coordinates of the lattice points of the array and δ is Dirac's delta function. In addition to the above differential equation we require that the integral of S 1 over the unit cell be zero. A Fourier series representation of S 1 and an efficient technique based on Ewald summation for evaluating S 1 are described by Hasimoto (1959). Substituting (2.15) into (2.14), and making use of (2.16), we find that the nondimensional pressure gradient is related to the sum of monopoles: where A 0 is the average monopole. The multipoles A α n andÃ α n and the constant U 0 in (2.15) are to be determined from the no-slip boundary condition u s = 0 on the surface of the tubes and (2.11), which states that the non-dimensional superficial velocity is unity. For this purpose it is convenient to re-expand u s around the centre of each tube. For example, u s is expanded near tube α as u s = ∞ n=0 u α n (r) cos nθ +ũ α n (r) sin nθ (2.18) with u α n (r) = a α n r −n + e α n r n (n > 1), The terms singular at r = 0 in the above expression arise from the singular part of S 1 at r = 0. Noting that S 1 behaves as −2 log r as r → 0 (Hasimoto 1959), and using the formulae for the derivatives of log r in Sangani & Yao (1988b), we obtain The coefficientsÃ α n are similarly related toã α n . The coefficients of the regular terms, such as e α n , are related to the derivatives of the regular part of u s at x = x α (Sangani & Yao 1988b). For example, where ξ n = n/4 for n > 2,ξ n = (n − 2)/4 for n > 3, and ξ 0 = ξ 1 =ξ 1 =ξ 2 = 0. In (2.21)-(2.22), u r s denotes the regular part of u s obtained by removing the singular part, −2 log r, from S 1 (x − x α ).
To determine the relation between U 0 in (2.15) and the multipole coefficients we must integrate u s over the area A s occupied by the shell-side fluid. Since the integrals of S 1 and its derivatives over the unit cell vanish, it is easier to evaluate the integral of u s over A s by integrating (2.15) over the unit cell and subtracting from it the integral of u s inside the tubes. With the non-dimensional superficial velocity taken as unity, the above procedure yields Care must be taken in carrying out the above integration to account for the singular nature of u α s at x = x α . Upon integrating, we obtain The last term on the right-hand side was missing in the expression given in Sangani & Yao (1988b). Fortunately, the omission of this term led to only a small numerical error in their results for the pressure drop. The no-slip boundary condition on the surface of the tube, together with the orthogonality of trigonometric functions, requires that Substituting for a α n and e α n from (2.20) and (2.21) into the expressions for u α n and applying (2.25) we obtain a set of linear equations in the multipole coefficients A α n . This set is truncated by retaining only the terms with n 6 N s to yield a total of 2N s + 1 equations in the same number of unknowns; solving it yields the velocity of the fluid on the shell side.

Determination of the concentration field
The solute concentration in the shell-side fluid is determined in a similar manner. A formal solution of (2.3) that is spatially periodic is given by where the spatially periodic function S 2 satisfies As shown by Hasimoto (1959) where the summation is over all reciprocal lattice vectors except k = 0. As mentioned earlier, Hasimoto (1959) describes a method for evaluating these functions using the Ewald summation technique. Substituting for C s and u s from (2.26) and (2.15) into (2.3) and using (2.16) and (2.27), we find that, in order for (2.26) to be the solution for C s , we must have To determine the multipoles B n , we use the same procedure as for determining the velocity field. Thus, we expand C s near the centre of each tube: and similar expressions forf α n . The coefficients appearing in the above expressions can be related to the coefficients A n and B n using the procedure outlined for the flow problem. Combining the above expansion for concentration with that inside the tube and applying boundary conditions of continuity of concentration and flux yields a set of linear equation for determining B α n , and hence the concentration field. The difference between the average solute concentration in the tube-and shell-side fluids is given by It is customary to write the overall mass transfer coefficient in terms of individual mass transfer coefficients on the shell and tube sides. The tube-side concentration drop is easily calculated to be where C w is the concentration at the surface of a tube. The tube-side Sherwood number is therefore This, of course, is the well-known result for the tube-side Sherwood number for the Graetz problem based on constant wall flux. We now define the shell-side Sherwood number Sh s via 1 For determining Sherwood numbers based on mixing-cup concentration differences (cf. (2.10)), we need to integrate the product u s C s over the area occupied by the shell-side fluid. This is difficult because it would require evaluating S 1 , S 2 , and their derivatives at many points outside the tubes. It is more efficient to solve instead for an auxiliary function ψ defined by Substituting for C s from (2.38) into (2.10) and using Green's theorem we obtain (2.39) The integral over ∂A s , which consists of the unit cell boundary and the surface of the tubes, vanishes owing to the boundary condition u s = ψ = 0 on the tube surface and the spatial periodicity of ψ and u s . On using (2.14) we obtain (2.40) A formal expression for ψ can be written in the same way as for u s and C s : (2.41) where S 1 , S 2 and S 3 , and their derivatives, are to be evaluated at x − x α , and ∇ 2 S 3 = S 2 . Expression (2.28) with m = 3 can be used to evaluate S 3 . The coefficients ψ 0 , E n and E n are to be evaluated from the boundary condition ψ = 0 on the surface of the tubes. Finally, since ∇ 2 S 1 = 4π/τ at all points outside the tubes, we require that The coefficients E α n can be determined in the same manner as the method used for determining A α n and B α n . The details may be found in Koo (2002). The mixing-cup-based concentration difference is expressed as (2.44) The expression for determining the shell-side concentration difference C s,c can be found in Koo (2002).

Results
Our results for the Sherwood number for periodic arrays with N = 1 were found to be in good agreement with those reported by Sparrow et al. (1961). We shall therefore present results for random arrays only. Figure 2 shows the inverse of the Sherwood numbers as functions of N, the number of tubes per unit cell, for φ = 0.1 and α c = ∞. The results were obtained by averaging the shell-side concentration difference over 100 hard-disk configurations for each N. A molecular dynamics code was used for generating hard-disk random configurations. We see that both Sh −1 s and Sh −1 s,c increase logarithmically with N. The solid lines in this figure indicate the slopes predicted by the theory, to be described next.
As mentioned in the Introduction the logarithmic divergence arises because the fundamental solution of the Laplace equation in a two-dimensional space is log r. The tube-side fluid acts as a source of solute while the shell side acts as a sink. Our theory will show that there is a net source due to the presence of each tube, and this implies that the concentration disturbance caused by a tube would grow logarithmically. To show this let us begin by deriving the equation for the conditionally averaged concentration, i.e. the ensemble-averaged solute concentration subject to a condition that a tube is present with its centre fixed at the origin. Inside the tube the solute concentration satisfies where C 1 and q 1 are, respectively, the conditionally averaged concentration and flux. Outside the tube, the equation governing the conditionally averaged concentration is slightly more complicated since a given point may lie inside another tube or outside all the tubes. Let χ be an indicator function whose value at a given point is unity if it lies inside a tube and zero otherwise. The conditionally averaged source density H 1 then is given by The apparent source due to the presence of a tube at the origin as seen from a distance R is therefore The above source must equal the net outward solute transfer from the surface r = R. At large r, the source density vanishes since the conditional averages converge to the unconditional averages and χu t 0 φ −1 = (1 − χ)u s 0 = 1. The integrand in the last integral in (3.3) therefore vanishes at large r and one may substitute R = ∞ for the purpose of evaluating the total apparent source due to the presence of a tube at the origin. Since C 1 (r|0) must be a function of r only for a random isotropic medium, and must satisfy the Laplace equation at large r where the source density H 1 vanishes, we must have that for large r Here, D * is the effective diffusivity at large r satisfying q 1 = −D * ∇ C 1 . The problem of predicting the effective diffusivity D * of a medium containing disks of diffusivity α c randomly dispersed in a medium of unit diffusivity has been examined by a number of investigators including Sangani & Yao (1988a) who presented results for D * as a function of α c and φ.
The behaviour of C 1 as predicted by (3.4) is valid for r large compared with unity (the tube radius) but small compared with the unit cell size, i.e. for 1 r h. On the unit cell length scale C 1 must, of course, satisfy the periodicity requirement. We analyse the problem using the method of matched asymptotic expansions with (3.4) representing the behaviour in the inner region, r h. In the outer region, valid for r = O(h), C 1 must satisfy Laplace equation to leading order, must be spatially periodic, and must match with (3.4) as r → 0. These conditions are satisfied by C 1 (r) = (Q ap /4πD * )S 1 (r) + const. (3.5) The constants in (3.4) and (3.5) need not be equal. It may be noted that the Laplacian of C 1 as given by (3.5) is not exactly zero since ∇ 2 S 1 = 4π/τ . Thus the apparent source at the centre of the tube is balanced by a uniform sink of strength Q ap /τ distributed throughout the unit cell. This sink strength, being O(h −2 ), does not affect the leading-order behaviour in the inner region. Now using the fact that S 1 → 2 log(h/r) + O(1) as r → 0, and noting that h 2 = πN/φ it is easy to show that (3.7) To determine the constant B we must evaluate the source density H 1 in (3.2) and integrate it over the space outside the fixed tube. The first term on right-hand side of (3.2) can be written as where P (r |0) = φg(r )/π is the probability density of finding a tube with its centre at r given the presence of a tube at the origin and g(r ) is the radial distribution function. The integration in (3.8) must be carried out over the area of a unit circle with |r − r | 6 1. The second term on the right-hand side equals the conditionally averaged velocity in a random array of fixed disks where u denotes the velocity field in a fixed bed of disks with u = u s for χ = 0 and u = 0 for χ = 1.
To determine the conditionally averaged velocity in a fixed bed of disks, we multiply the momentum equation (2.14) by the fluid indicator function 1 − χ and ensemble average the resulting expression with a disk fixed at the origin: (3.10) Inside the disks the velocity is zero and hence Taking the Laplacian we obtain Note that the second term on the right-hand side of the first equality vanishes owing to the no-slip boundary condition. Since ∇χ equals the unit normal vector pointing into the disk multiplied by a delta function centred at the disk circumference, the last term on the extreme right-hand side of (3.12) is given by − ∇χ · ∇u 1 (r|0) = |r−r |=1 P (r |0)n · ∇ u 2 (r|0, r ) dr .
( 3.13) The integral can be expressed in terms of an integral on a particle centred at r by use of a Taylor series expansion to yield − ∇χ · ∇u 1 = P (r|0) (r − r)n · ∇ u 2 (r |r, 0) dr + . . . (3.14) The first term on the right-hand side can be evaluated from the local expansion of the velocity field near a representative tube α (cf. (2.18)) which shows the integral to equal πG + 4πA α 0 . The second term is related to the stresslet induced by the presence of the disk and contributes to the effective viscosity of the medium. Ignoring the higher-order terms in (3.14), and combining (3.10) and (3.12), we obtain where µ B is the Brinkman viscosity and g(r) is the radial distribution function defined by P (r|0) = ng(r), n being the number density of tubes. The Brinkman viscosity is defined via the closure µ B (r)∇ u 1 = ∇ u 1 (r|0) + P (r|0) |r−r |=1 (r − r)n·∇ u 2 (r |r, 0) dr .
3.1. Effective-medium approximations 3.1.1. Velocity field To make further progress in determining the conditionally averaged velocity we employ an effective-medium approximation in which a simple, single disk model is used. Since the disks are non-overlapping, there is a net depletion of the number density of disks in the immediate vicinity of a disk. To account for this we assume an exclusion region around the fixed disk at the origin for 1 < r < R. The conditionally averaged velocity satisfies the equations of motion for a single-phase flow (cf. (2.14)) in this exclusion region. Outside this region the fluid-disk medium is replaced by a medium whose properties are consistent with the average properties and average equations of motion for flow through a fixed bed of disks. The radius R is chosen so that the exclusion area around a fixed disk equals that in the actual bed. Since the reduction in number density near the disk equals n − P (r|0), we require that The quantity on the right-hand side can be expressed in terms of the zero wavenumber limit of the structure factor S(0) to yield where the use has been made of the relation φ = nπ. The above choice for the exclusion radius R was made by Dodd et al. (1995) who showed that the effectivemedium approximation based on this value of R agrees very well with the results of rigorous computations for the mobility of integral membrane proteins in bilipid membranes, modelled as suspensions of disks. Subsequently, Wang & Sangani (1997), Sangani & Mo (1997), and Spelt et al. (2001) have used a similar model for estimating properties of suspensions of disks and spheres. The computational results presented in figure 2 corresponded to hard-disk random configurations. S(0) for these configuarions can be evaluated using (Chae, Ree & Ree 1969) S(0) = (1 − 1.9682φ + 0.9716φ 2 ) 2 1 + 0.0636φ − 0.5446φ 2 − 0.4632φ 3 − 0.1060φ 4 + 0.0087φ 5 . (3.19) Note that S(0) → 1 − 4φ as φ → 0 and this yields the exclusion radius R equal to 2 in the limiting case of very dilute random suspensions. Next, we introduce the following closure relation for the monopole in the effective medium: 4πn A 0 2 (r) = −κ 2 u 1 (r).
(3.20) Finally, we take µ B = 1 for r < R and a constant equal to µ * for r > R. We also take g(r) = 0 for r < R and g(r) = 1 for r > R. The conditionally averaged velocity therefore satisfies ∇ 2 u 1 = −κ 2 u 0 for 1 < r < R, (3.21) µ * ∇ 2 u 1 = κ 2 ( u 1 − u 0 ) for r > R. (3.22) The above equations together with the boundary conditions u 1 = 0 at r = 1, u 1 → u 0 = 1 as r → ∞, and the continuity of u 1 and tangential stress at r = R complete the description of the effective-medium model for determining the conditionally averaged velocity. The above equations can be solved readily to yield with κ * = κ/ √ µ * . Here, K 0 is the modified Bessel function of zeroth order. The constants β and A 0 are to be determined using the continuity of velocity and stress at r = R. The Brinkman viscosity µ * is related to the effective diffusivity in reacting media. Calculations of Sangani & Behl (1989) for diffusion into a semi-infinite medium of reacting spherical traps seem to indicate that the diffusivity of reacting media is close to the diffusivity of the medium surrounding the traps. This is equivalent to choosing µ * = 1.
To test the validity of the conditionally averaged velocity field obtained from the above model, we computed directly the conditionally averaged velocity field in a fixed bed of disks from the results of our rigorous numerical simulations. The velocity was evaluated at all the points of a 10 × 10 grid in the unit cell. For a simulation with N disks per unit cell, this provides a total of 10 2 N data points for the conditionally averaged velocity versus the distance r from the centre of a disk. Simulations were carried out for 100 configurations with N = 100 to yield a total of 10 6 data points. The result for φ = 0.1 is shown in figure 3. In these simulations the unit cell size h is approximately equal to 56. Since g(r) would be close to unity for, say, r > 4, we can expect (3.23) to be reasonably accurate for r > 4. Regression analysis of the data for 4 < r < 7 gave µ * = 1.063 while that for 3 < r < 8 gave µ * = 0.951. We see that in both cases the Brinkman viscosity is not very different from unity. Figure 3 shows a comparison between the results of exact computations of the conditionally averaged velocity and the velocity predicted by the effective-medium model for several selected values of µ * . We note that the velocity profile is relatively insensitive to µ * and that the predicted velocity profile agrees very well with the computed one.
The closure (3.20) for A 0 2 can also be verified using the results of numerical simulations. A simulation with N disks has N(N − 1)/2 pairs of disks and give that many data for the conditionally averaged monopole A 0 2 . Assuming that the closure (3.20) applies, we can estimate the conditionally averaged velocity from the data for A 0 2 . The results obtained with N = 100, and 100 configurations with φ = 0.1 are also shown in figure 3. We see good agreement for r greater than about 3. Equation (3.18) gives R = 1.84 for φ = 0.1. The permeability predicted using the effective-medium approximation with this value of R equals 5.16 while the simulations gave 5.34, about 4% greater. Results for other values of φ are shown in figure 4. We see that agreement is good even for high values of φ. The effective-medium curve shown in figure 4 was obtained by taking the Brinkman viscosity to be unity, i.e. µ * = 1.

Concentration field
We now develop an approximate theory for the conditionally averaged solute concentration field. Using the closure relation that relates the solute flux to the concentration gradient, (3.2) reduces to where g(r) = P (r|0)/n is the radial distribution function, n being the number density of tubes. As in the case of the approximation for the conditionally-averaged velocity, we take D = 1 and χ = 0 for 1 < r < R and D = D * and χ = 1 for r > R, where, as mentioned earlier, D * is the effective diffusivity of a medium consisting of disks of diffusivity α c suspended in a medium of unit diffusivity. Numerical results for D * as a function of α c and φ, the area fraction of the disks, have been reported by a number of investigators. We shall use the results for random arrays of disks presented by Sangani & Yao (1988a). On solving (3.24) we obtain C 1 = κ 2 (r 2 /4 − r 4 /16)/4 + Ar 2 (1 − log r)/2 + E * log r + H * for 1 6 r 6 R βK 0 (κr)/(D * κ 2 ) + E log r + H for r > R.
(3.25) The expression for r > R applies to distances that are small compared with the unit cell size h. For distances comparable to h, an outer region approximation for C 1 can be obtained by noting that, since u 1 u 0 = 1 for r = O(h), C 1 satisfies the Laplace equation. It is easy to show therefore that (3.26) The unconditionally averaged solute concentration C 0 can be set to zero with no loss of generality. The relation among B, E and H can be determined by requiring that (3.26) agrees in the limit r/ h → 0 with that for the inner region given by (3.25) as r → ∞. For small r/ h, S 1 for a square lattice is given by Matching the solution for C 1 in the two regions therefore yields (3.28) The constants B, E * , H * and C t can be determined now from the continuity of concentration and flux at r = 1 and r = R. Since we have taken C 0 = 0, the concentration difference C = C t − C s can be evaluated using C s = −φ C t /(1 − φ). Figure 5 shows a comparison between the conditionally-averaged concentration determined by the above approximate theory and the rigorous results obtained from numerical simulations with N = 100, α c = ∞, and φ = 0.1. The rigorous results indicated by filled circles were obtained by computing C on a 10 × 10 grid placed in the unit cell as in the case of the conditionally-averaged velocity. Alternatively, the conditionally-averaged concentration can also be estimated from the conditionally-averaged concentration in the tubes, C t 1 (r|0) using the approximation C 1 (r|0) = (1 − φ) C + C t 1 (r|0) if we assume that C s 1 − C t 1 is independent of r. The results obtained in this way are shown by open circles. We see a good agreement between the two, indicating the validity of the assumption that the difference between the conditionally-averaged solute concentration and the conditionally-averaged concentration in the tube is independent of r. The thin line in the figure represents the inner-region effective-medium approximation while the thick line represents the outer region. The outer-region solution depends on r and θ, the angle between the vector r and the x 1 -axis. The results shown in figures were obtained by averaging over θ. We see that effective-medium theory gives a fairly accurate prediction of the conditionally-averaged concentration.
The results of the effective-medium approximation can be expressed as with λ 1 = 2φB/(1 − φ). Figures 6 and 7 present the results for λ 1 and λ 2 as functions of φ. The filled circles in these figures are the estimates of λ 1 and λ 2 obtained by fitting the results of numerical computations of Sherwood numbers as functions of N. We see that λ 1 is reasonably well predicted. Although the agreement for the O(1) coefficient, λ 2 , does not appear as good, it must be noted that variations and the magnitude of λ 2 are relatively small. Indeed, as seen in figure 8, which compares the simulation results for the inverse Sherwood number as a function of N with the predictions of the effective-medium theory, the agreement with the theory is reasonably good. Figure 9 shows λ 1 as a function of α c , the ratio of tube-to shell-side solute diffusivities. We see that λ 1 varies significantly with α c , The apparent source of solute Q app is approximately the same but D * increases significantly as α c is increased, and this in turn causes λ 1 to decrease as α c is increased.

Product of velocity and concentration fields
An approximate theory for determining the conditionally-averaged ψ, and hence the Sherwood number based on the mixing-cup solute concentrations in the shelland tube-side fluids, can be developed in a similar manner. Thus, it can be shown that the conditionally-averaged ψ satisfies (Koo 2002)  As in the case of the conditionally-averaged concentration, this equation is satisfied in the inner region, r h. In the outer region, ψ can be shown to be given by (3.31) 10 -2 10 -1 10 0 10 1 10 2 10 3 λ 1 Figure 9. λ 1 as a function of α c for φ = 0.4.
(3.32) The unknowns ψ 0 and γ are determined by requiring that both ψ 1 and its derivative are continuous at r = R. Finally, the mixing-cup solute concentration in the shell-side fluid is determined using C s,c = G ψ 0 = 4πnκ 2 ψ 0 (cf. (2.40)). Figure 10 compares the prediction of the above theory with the simulation results for ψ 1 obtained by determining ψ on a 10 × 10 grid placed inside a unit cell with φ = 0.1. The figure also shows the results for the conditionally-averaged ψ obtained by combining the data for monopole E 0 2 as a function of r with the closure 4πn E 0 2 = κ 2 ( ψ 0 − ψ 1 (r|0) ( 3 . 33) used in deriving (3.30). We see that the theory is in very good agreement with the simulation results, especially at larger values of r. For small r, however, the theory and simulation results show considerable discrepancy. The discrete nature of the disksfluid system is apparently not captured well by the effective-medium approximation for small r. Figure 11 compares the numerical simulation results for the inverse Sherwood number based on mixing-cup concentrations with that predicted by the theory for large N. The slopes of the lines shown in the figure depend on the coefficient of the log N term, and are seen to be in reasonably good agreement with the simulation results. The O(1) constant on the other hand is underpredicted by the theory. This may be due to a discrepancy between the theory and simulation results for ψ 1 − ψ 0 at small r.

Effect of the shell boundary
We have analysed so far the problems of determining the concentration and velocity fields in an unbounded periodic array consisting of N tubes per unit cell. The case of greater practical significance is that of N tubes surrounded by a shell boundary. Since the conditionally-averaged velocity disturbance approaches the unconditionally averaged velocity within a distance comparable to the tube radius, due to Brinkman screening, the effect of the shell boundary will be insignificant as long as the shell dimensions are much greater than the tube radius. However, the conditionallyaveraged concentration, being divergent on the lengthscale of tube radius, will be affected by the shell boundary. The results obtained for the periodic arrays must therefore be modified to account for the shape of the shell boundary. The coefficient of the leading, O(log N), term in the inverse Sherwood number, i.e. λ 1 in (3.29), will be the same as in the case of a periodic array but the O(1) coefficient, λ 2 , will be affected by the shape of the boundary. We shall outline here the problem that must be solved to determine the influence of the shell boundary shape and give a detailed analysis for the special case of a circular shell of radius S (cf. figure 1).
The inner-region analysis for the conditionally-averaged concentration remains the same as before but the outer-region analysis must be modified. In lieu of (3.26) we now have where we have replaced the Green's function for Laplace's equation in periodic arrays, S 1 , by a Green's function G defined below. It should be noted that, unlike the case of periodic arrays, the conditionally-averaged concentration is now also a function of the position vector r of the centre of the tube. The constant B, being related to the net source due to the presence of the tube, is the same as before. The Green's function is defined via where τ is the cross-sectional area of the shell, ∂τ is the boundary of this cross-section, and n is the unit outward normal at the boundary. The matching of the inner-and outer-region expressions for the conditionallyaveraged concentration requires taking the limit r → r of (3.34). Let G r (r ) = lim r→r [G + 2 log |r − r |]. (3.38) Note that G r is simply the regular part of G at r = r . Now the matching process is the same as in the case of periodic arrays and the resulting relations among E, B and H are obtained by replacing S r 1 ≡ 2 log h − 2.6232 in (3.28) by G r . The difference in the tube-and shell-side fluid concentrations is therefore directly related to G r . To determine the overall shell-side mass transfer coefficient we must average G r over the entire cross-section τ . Let The Sherwood number for the shell-and-tube configuration (st) is then related to that for periodic arrays (per) by An analytical solution for G is possible for the special case of a circular shell. It can be shown that the solution of (3.35)-(3.37) is given by where S is the shell radius non-dimensionalized by the tube radius a and θ and θ are the angles between vectors r and r and the x 1 -axis. In deriving the above result use has been made of the following expansion for r > r : and a similar identity obtained by interchanging r and r for r > r. The regular part of G at r = r is given by (3.43) On averaging this over the shell cross-section we obtain G r = 2 log S + 3/2.
(3.44) Now using S 2 = N/φ = h 2 /π in (3.40) we obtain (3.45) indicating that the Sherwood number for the shell-and-tube configuration is lower than that for the peridic array having the same N. The concentration difference between the tube-side and the surrounding shell-side fluid is greater for the tubes that are close to the shell boundary than in the centre because of the no-flux condition at the shell boundary and as a result the overall Sherwood number is smaller.

Mass transfer at small axial distances
The analyses in the previous sections showed that the mass transfer coefficients for fully developed laminar flow conditions are very small. From the practical point of view then the contactors must be made short enough so that the concentration boundary layers around the tube surface remain thin compared with the average spacing between the tubes. This condition is approximately satisfied when L/a Pe 1/3 , L being the length of the tubes. In this section we first determine the mass transfer coefficient for short exchangers satisfying the above condition and then propose a correction factor to account for the overlap of concentration boundary layers in random arrays.
We shall limit our analysis to the case when the concentration of the tube wall is constant for x 3 > 0. It is relatively straightforward to show that the shell-side Sherwood number for this case can be evaluated using where τ w is the non-dimensional traction at the surface of a tube. Note that since the concentration boundary layers are very thin, the Sherwood number based on mixing-cup concentration difference is the same as that based on area average. An estimate of the traction at the tube surface may be obtained either using the cell theory for periodic arrays or using the effective-medium theory for random arrays. In both cases the velocity field is independent of θ and the wall traction is related to the monopole A 0 , which, in turn, can be related to the permeability of the array. The final result is where K is the permeability of the array non-dimensionalized by a 2 . Table 1 compares the predictions obtained from the above theory with the results of numerical simulations in which the wall traction on each tube was calculated and the integral in (4.1) evaluated using the Simpson's rule. We see a good agreement between the two. The results of numerical simulations are also seen to be in good agreement with those presented by Bao et al. (1999) who considered both the constant wall concentration and flux cases. For random arrays the assumption of non-overlapping concentration boundary layers is too restrictive and the above result must be modified before it can be applied to practical situations. Figure 12 shows the results of numerical simulations for determining F (δ, φ), the average of the fraction of a tube surface that is within a distance of 2aδ from the surface of another tube. The simulations were obtained by averaging over 100 configurations with N = 100. As one would expect F → 1 as δ → 0 and F decreases more rapidly with δ as φ is increased. Thus, the probability of concentration boundary layers overlapping in finite-length contactors increases with increasing φ.
The solid lines in figure 12 correspond to fitting the results of numerical simulations according to (4. 3) The results for δ 0 for 0.1 6 φ 6 0.5 can be fitted according to (4.4) Since the regions in which the concentration boundary layers begin to overlap do not contribute to the mass transfer as effectively as the regions in which the layers do not overlap, we propose that the formula (4.2) be corrected by multiplying its right-hand side by F with aδ representing the average concentration boundary layer thickness.
The concentration boundary layer thickness, defined as the distance from the surface of a tube over which 95% of the concentration drop occurs, is given by δ = 2.37 x 3 τ w P e 1/3 . (4.5) This thickness is, of course, a function of angular position on the tube surface since the wall stress is not uniform. We replace τ w by its average value which can be estimated from (4.1) as τ 1/3 w = β z (φ)/0.5384. Thus, we propose that δ to be used in the correction factor F be taken as (4.6) with A z = 0.5384 × 2.37 = 1.28. Note that in terms of this definition of δ, the Sherwood number can be alternatively expressed as This expression is particularly convenient for determining the average Sherwood number for a contactor of length aL. The particular form chosen for fitting the results of numerical simulations for F , i.e. (4.3), allows rather straightforward integration of the above expression for the local Sherwood number to yield the average Sherwood number: This can be added to the Sherwood number corresponding to the fully developed case to yield an estimate for arbitrary values of L.

Comparison with experiments
The most commonly cited correlations for Sh in hollow-fibre contactors are due to Yang & Cussler (1986), Prasad & Sirkar (1988), and Costello et al. (1993). Yang & Cussler (1986) carried out experiments for φ = 0.03, 0.26, and 0.4 but summarized the results in their table 2 for only φ = 0.03 and 0.4. The latter corresponded to a module of 2100 fibres with a = 0.022 cm and aL = 22 cm enclosed in a rectangle-shaped shell. The non-dimensional length L was therefore 1000. The experiments were carried out with water containing dissolved oxygen on the shell side and either vacuum or nitrogen gas on the tube side. The fibres were made of hydrophobic membranes which offered very little resistance to the mass transfer of oxygen as the pores of the membrane remained relatively dry permitting easy diffusion of oxygen. The Schmidt number for the oxygen-water system is 476 and, since the diffusivities in the gas phase are much greater than in the liquid phase, we may take α c = ∞. The experiments were carried out at several flow rates with the Reynolds number in all cases less than 1 and the Sherwood number was found to be independent of the flow rate and equal to 0.08. (Note that the value 0.24 reported in their table 2 corresponds to using the hydraulic diameter d h = 2a(1 − φ)/φ as the characteristic length, which for φ = 0.4 equals 3a.) The investigators were apparantly surprised by their result as the following quotation from their paper suggests. "Under many circumstances we would interpret this as evidence that the membrane resistance has become important. In this case, we are not so sure. The membrane resistance almost certainly is not important for less densely packed fibres, so it seems strange that it is suddenly significant here. Instead, we suspect that there is major channelling through the closely packed fibres, and that the mass transfer is controlled by diffusion through nearly stagnant liquid trapped between the fibres, not by the membrane".
With Re less than 1 and Sc = 476, the Péclet number P e is less than 476, which is small compared with L of 1000 in the above experiments, and we expect the results for the fully developed case to be applicable. The Sherwood numbers for N up to 100 were given in figure 11 for the case α c = ∞ corresponding to tube-side diffusivity much greater than the shell side. Extrapolation to N = 2100 and φ = 0.4 yields Sh of about 0.09, quite close to the value reported by Yang & Cussler considering that the rectangular shaped shell used in their experiment is not the same as the periodic unit cell used in our computations. Note that these values for random arrays are much smaller than a Sherwood number of about 4 for a square array with φ = 0.4. This comparison also illustrates the significance of the shell-side mass transfer resistance in these contactors when a liquid is flowing on the shell side and a gas on the tube side. Prasad & Sirkar (1988) carried out experiments with liquids flowing on both tube and shell sides, the liquid on one side being typically water and on the other side an organic liquid (xylene, n-Butanol, or methyl isobutyle ketone). The mass transfer rates of a solute (acetic acid, succinic acid, or phenol) across both hydrophobic and hydrophilic membranes were determined for various φ, Re and L. Their results apparently can all be collapsed onto a single generalized chart shown in their figure 6 to give the correlation Re 0.66 Sc 0.33 n (5.1) with A = 11.7 and n = 0.78. (Note that the values reported by these investigators, i.e. A = 5.7 and n = 1, and quoted by several subsequent investigators do not quite correspond to the solid line shown in their figure.) They used both hydrophilic and hydrophobic membranes and found that the shell-side Sherwood numbers obtained using a hydrophilic membrane were slightly greater than those using a hydrophobic membrane although the scatter in data in both cases is significant. Detailed comparison with their data is difficult because the results obtained with different Sc, a/L and φ were all plotted in the same figure without indicating which data were obtained with which system. They suggest that the correlation (5.1) is valid for Re < 250φ/(1−φ) (less than 500 based on hydraulic diameter) and 300 < Sc < 1000. Figure 13 compares the experimental correlation with the theory prediction for φ = 0.4 and 0.1 with Sc = 350 and L/a = 1100, corresponding approximately to the Xylenewater-acetic acid system. The experimental correlation for these two area fractions is indicated by diamonds and filled circles. These data correspond to Reynolds numbers based on hydraulic diameter of less than 500. The vertical bars correspond roughly to the scatter in the data shown by Prasad & Sirkar. The solid lines correspond to the theoretical prediction corresponding to choosing δ as roughly the distance from the tube wall where the concentration drops by 95% of the total drop. We see a very good agreement between the theory and experimental correlation. The dashed lines in this figure correspond to taking the correction factor F to be unity, i.e. not accounting for the overlap of the concentration boundary layers. We see that it overestimates the mass transfer coefficients considerably for the range of parameters considered in the experiments. It may be noted that the actual values of the Sherwood numbers in these experiments are rather small, O(1), and much lower than would be expected for periodic arrays.
Prasad & Sirkar also plotted the data of Yang & Cussler for φ = 0.26 and showed that their data were roughly in agreement with the correlation given by (5.1) with perhaps slightly greater n when the Reynolds number based on hydraulic diameter is less than 500. These data agree well with the theory also. Yang & Cussler presented data to much higher Reynolds numbers and these higher Reynolds number data deviate strongly from Prasad & Sirkar's correlation and our theory. This is also true of their data for φ = 0.03 obtained with a bundle of 16 fibres. Sherwood numbers obtained with φ = 0.03 are rather large, O(10), and it is possible that the assumption of fully developed laminar flow may not be appropriate for this system. Large spacing between the fibres would require longer distances for the flow to be fully developed. The Sherwood number for developing flow is expected to be higher as thin hydrodynamic boundary layers in the developing flows contribute to higher mass transfer rates than in fully developed laminar flows. Also, the end effects of introducing the shell-side fluid perpendicularly to the fibre bundle may persist over longer lengths. Finally, Costello et al. (1993) obtained Sherwood numbers for φ, a/L and Re values that are comparable with the values considered by Prasad & Sirkar (1988) but obtained much higher values of Sherwood numbers. They also measured the shell-side pressure drop and found it to increase more rapidly with the velocity than predicted by laminar flow. It appears that the longitudinal, laminar flow assumption is not appropriate for their system.

Summary
We have determined mass transfer coefficients for longitudinal flow through a periodic array of N randomly placed tubes. The concentration disturbance caused by each tube grows logarithmically with the distance for distances small compared with the unit cell size. This leads to large concentration differences between the tubeand shell-side fluids for systems with a large number of tubes per unit cell (or large shell to tube radius ratio), and, consequently, small Sherwood numbers. Effective-medium theory appears to not only give reasonable estimates of the overall transport properties such as the Sherwood numbers and the permeability, but also provide reasonably accurate conditionally-averaged fields.
Financial support for this work was provided by the National Science Foundation under Grant No. CTS-9909234). The computations were performed using the resources by the National Center for Supercomputing Applications at University of Illinois at Urbana-Champaign.