- Email: [email protected]
Techniques for computational fluid dynamics modelling of flow in membrane channels
Journal of Membrane Science 211 (2003) 127–137
Techniques for computational fluid dynamics modelling of flow in membrane channels Dianne E. Wiley a,∗ , David F. Fletcher b a
UNESCO Centre for Membrane Science and Technology, School of Chemical Engineering and Industrial Chemistry, University of New South Wales, Sydney 2052, NSW, Australia b Department of Chemical Engineering, University of Sydney, Sydney 2006, NSW, Australia Received 30 December 2001; received in revised form 15 May 2002; accepted 5 September 2002
Abstract Accurate modelling of the flow and concentration polarisation in pressure driven membrane processes is inhibited by the complex couplings in the flow equations along with any added effects of variable solution properties. A generic computational fluid dynamics (CFD) model has been developed which incorporates these effects and describes the flow across the membrane wall. The results have been validated against classical solutions available in the literature. Extended work indicates that overly simplified expressions for the dependence of viscosity and diffusivity on concentration produce velocity and concentration profiles that may grossly misrepresent reality. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Concentration polarisation; CFD; Theory; Ultrafiltration; Reverse osmosis
1. Introduction The modelling of flow and concentration polarisation in membrane or thin channels with permeable walls is not new. Numerous approximate one-dimensional [1–6] and two-dimensional [7–38] models of varying degrees of complexity and simplification have appeared in the literature. A succinct review of much of the early work was prepared by Kleinstreuer and Belfort [39]. The attraction of these analytical and semi-analytical models is the ability to tailor them for relatively quick and easy investigations of the effect of different approximations and parameters on the performance of membrane ∗ Corresponding author. Tel.: +61-2-98385-4304; fax: +61-2-98385-5966. E-mail address: [email protected] (D.E. Wiley).
systems. Models exist which permit investigation of the effects of variable permeation [17–23], changes in rejection [17,18], effects of wall slip [24,25], variable physical properties of the solution being processed [26–36] and gravitational effects [37,38]. Despite the simplifications, many of the models have proved to be remarkably accurate at predicting membrane performance under some limited circumstances. Development of a generic model with widespread applicability to membrane systems requires much more rigorous and robust techniques than those described above. Computational fluid dynamics (CFD) provides one such technique. There have been many fewer attempts at modelling membrane systems in this manner. This is most likely due to the complexity of the problem and the fact that the ODE based approaches described previously are computationally very efficient. However, CFD models are needed if
0376-7388/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 6 - 7 3 8 8 ( 0 2 ) 0 0 4 1 2 - X
128
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Nomenclature B0 B1 B2 c C+ D h H L p R Re Sc u v x X y
constant in Eq. (8) constant in Eq. (5) constant in Eq. (5) concentration (wt.%) dimensionless concentration, c/c0 diffusivity (m2 s−1 ) half flow channel height (m) half permeate channel height (m) channel length (m) pressure (Pa) rejection coefficient Reynolds number Schmidt number velocity in the axial direction (m s−1 ) velocity in the transverse direction (m s−1 ) axial direction (m) fractional water removal, vw0 x/u0 h transverse direction (m)
Greek letters α D/vw h β constant in Eq. (4) µ dynamic viscosity (Pa s) ν kinematic viscosity (m2 s−1 ) ρ density (kg m−3 ) ξ 3B13 x/ h Sc Re Subscripts p permeate w wall 0 inlet
the flow in complex geometries, which include recirculation, are to be studied. A variety of techniques have been used and are summarised below. Huang and Morrissey [40] used the finite element technique to solve for the concentration field given an approximate solution to the velocity field, by neglecting stream-wise diffusion and assuming constant properties. This work is very limited but it highlighted the need for extremely high resolution near the membrane surface because of the extremely thin concentration boundary layer. Lee and Clark [41] performed a similar analysis but this time solved the concentration equation using a finite difference method. They also
included a shear induced diffusion term, as well as the usual molecular diffusion term, in the concentration equation. In addition, they developed an expression for the permeate flux which is based on diffusion through a stagnant cake layer assumed to develop on the membrane surface. Both works are limited by the fact that they assume the velocity field to be known. Hansen et al. [42] used a spectral method to solve the laminar Navier-Stokes equations for channel flow with permeate flow through a wall. They allowed the permeate velocity to be a function of the static and osmotic pressure differences and made the viscosity depend on concentration. The spectral method used had high spatial accuracy and required relatively few elements but used polynomials of degree five for concentration and pressure and three for velocity, giving a large number of degrees of freedom. Whilst the method worked well for the parabolic flow considered and for the case of 100% rejection, it is not clear how easily it can be extended to complex geometries and to cases where the rejection is less than 100%. Geraldes et al. [43,44] performed finite volume simulations for permeation of water through a membrane for a simple channel flow. They solved both the flow and the concentration fields for a constant property situation. The simulations used a solute of salt, so that the diffusivity was relatively large. Despite this they needed very fine meshes to obtain a grid independent solution. They present considerable discussion on differencing schemes but their conclusions are of limited value, as they apply to a flow-aligned geometry with a relatively thick concentration boundary layer. Later work [45] has extended this model to include variable fluid properties. The model has been validated against experimental data on nanofiltration and used to fit effective rejection coefficients. There are almost no simulations of flow in thin channels or membrane systems in three dimensions. Recent work by Karode and Kumar [46] using relatively coarse grids and first order differencing for convection has investigated the effect of spacers on pressure drop and shear stress, but as yet does not include the effect of permeation on flow and concentration polarisation. The present work is directed towards removing some of the assumptions discussed above and to validate a general purpose finite volume model for use in the simulation of membrane systems. As a first step, we performed simulations in two dimensions
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
to enable validation of the modelling framework but note that this approach is easily extended to certain three-dimensional flows.
2. Theory The starting point for the work presented here is the commercially-available finite volume code CFX4 (CFX, 1998) which solves the conservation equations on a co-located grid using the SIMPLEC algorithm [47]. Specifically, we have assumed the flow to be steady, incompressible with constant density, and laminar. The transport properties µ (viscosity) and D (diffusivity) can be a function of the concentration. Given these assumptions, the equations for conservation of mass, momentum and solute concentration are: ∂u ∂v + =0 ∂x ∂y
(1)
∂u ∂u ρ u +v ∂x ∂y ∂p ∂ ∂u ∂ ∂u =− + µ + µ ∂x ∂x ∂x ∂y ∂y ∂v ∂v ρ u +v ∂x ∂y ∂p ∂ ∂v ∂ ∂v =− + µ + µ ∂y ∂x ∂x ∂y ∂y
When the concentration polarisation is large, it is important that this condition is applied well downstream of the end of the membrane, otherwise the solution over the last part of the membrane can be affected significantly. At the symmetry plane, v = 0, and the normal gradients of the tangential velocity (u) and the concentration (c) are set to zero. At the walls, the no slip condition is imposed and the normal derivative of the concentration is set to zero. At a membrane, the conditions are more complex, as flow permeates through the wall. For a membrane located at the plane y = 0, the tangential velocity (u) is set to zero and the normal velocity is set to vw . In some cases the wall (permeate) velocity was set to a constant value and in others it was made a function of the wall concentration of solute. The expression used was either that of Brian [17], which takes the form cw vw = v0 1 − β (4) c0 or that of Doshi et al. [30] in which v0 B1 cw 1 − B2 vw = h Sc c0
(2)
and ∂c ∂ ∂c ∂ ∂c ∂c = D + D (3) u +v ∂x ∂y ∂x ∂x ∂y ∂y 2.1. Boundary conditions
129
(5)
The boundary condition on the concentration results from a balance of the convective and diffusive fluxes. In addition, it must take into account the fact that not all of the solute permeates through the membrane. This is done via the use of a rejection coefficient, R, and the concentration boundary condition is given by D
∂c + vw Rcw = 0 ∂y
on the feed side of the membrane and ∂c D + vw cp = (1 − R)cw vw ∂y
(6)
(7)
on the permeate side of the membrane. At the inlet, the flow is assumed to be fully developed and a parabolic flow is specified. As this boundary condition is applied at the start of the membrane, then a transverse velocity, as well as a normal velocity is specified, based on the solution of Berman [7]. Other arbitrary inlet flow conditions could be specified for different applications with appropriate modification of the inlet channel length if required. A uniform inflow concentration of c0 is specified. At the outlet, all derivatives in the flow direction are set to zero.
2.2. The geometry and computational mesh used A model, shown schematically in Fig. 1, was set up to represent a section of a channel with a membrane separating two regions. The physical dimensions of the geometry were varied to match the conditions in the cited references. For example, in the case of the comparison with the data of Brian [17] the channel was 2.7 m long, the feed channel had a half channel
130
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Fig. 1. Geometrical model of the flow channel used for the simulations.
height of 1 mm and the permeate channel had a height of 0.6 mm. Flow is imposed in the top section of the channel and occurs in the bottom channel as mass is transferred through the membrane. Typically a mesh with 50 cells across the half channel, 10 cells across the permeate channel and 60 cells along the channel were used. A non-uniform mesh was used, with the mesh density being higher at the start of the membrane and near to the surface of the membrane. Care was taken to avoid excessively high cell stretch rates. In all cases, we checked the grid independence of the solution and made sure that there were sufficient cells in the concentration polarisation layer that the thickness of this layer and the peak concentration were independent of the mesh used. 2.3. Numerical procedure All convective terms were differenced using a MUSCL scheme with the van Leer limiter to provide high spatial accuracy without the introduction of unphysical dispersive errors [48]. Standard second order differencing was used for the diffusive terms. Stone’s method was used to solve the resulting matrix equations for the velocity and concentration fields and incomplete Cholesky conjugate gradient (ICCG) was used for the pressure. Under-relaxation factors of 0.65 on the velocity components and unity on the pressure and concentration worked well for situations in which the polarisation was not great, such as those simulated by Brian [17]. However, when very high concentrations developed near the wall, resulting in a very large fluid viscosity, it was necessary to reduce the relaxation factors on all but the pressure equation to account for the extreme non-linear coupling
between the viscosity and the solute concentration. In extreme cases, values as low as 0.05 were required, resulting in the need to perform several thousand iterations to obtain convergence. In some circumstances we also found it to be beneficial to under-relax the viscosity. However, in all cases we were able to obtain extremely good convergence, with the total mass residual being reduced to 0.05% of the inflow mass.
3. Results and discussion The model developed in this paper has been validated against a number of ‘classical’ semi-analytical solutions that are well accepted in the membrane literature. Comparisons to be shown encompass the effects of various combinations of variable and constant wall flux, variable and constant solution properties and several values of constant rejection. Unfortunately there is no single solution available in the literature that covers all of the combinations. A range of validations of increasing complexity has therefore been selected to demonstrate the general validity and utility of the model developed here. In comparing the data, it must be remembered that the CFD model performs a rigorous evaluation of the complete transport equations presented in Eqs. (1)–(3) while the semi-analytical methods each approximate the solution of the equations in different ways. Thus, some discrepancies between the various methods would be anticipated. This will be elaborated on in the ensuing discussion. In order to demonstrate that the boundary conditions at the membrane have been correctly specified, Fig. 2 compares our CFD solutions with data published by Brian [17] for variable wall flux (α = 0.27, β = 1),
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
131
Fig. 2. Comparison of the CFD data with the data of Brian for variable wall flux (α = 0.27, β = 1), constant solution properties and constant rejection (100 or 90%).
constant physical properties of the solution (ρ, D, µ) and constant rejections of 100 or 90%. Brian [17] used a finite difference technique with truncated versions of the velocity profiles determined by Berman [7] using a perturbation technique. Both Brian [17] and Berman [7] demonstrated the acceptability of truncating the velocity profiles when only small amounts of fluid were removed over the length of the permeating channel. Fig. 2 shows the excellent numerical accuracy of the model developed, although, it consistently predicts slightly smaller values of the wall concentration. The CFD model solves the full Navier-Stokes and concentration equations, whereas Brian [17] uses a truncated velocity distribution and ignores stream-wise diffusion of the solute. Comparison of the calculated velocity profile with the truncated expression shows a difference of less than three percent in the wall region for the cases investigated in this paper. This difference, together with the neglect of stream-wise diffusion, appears sufficient to explain the minor differences. Fig. 3 compares the effect of some approximate analytical solutions presented by Doshi et al. [30] against the rigorous CFD solution for constant wall flux (B2 = 0), constant physical properties (B0 = 0) and 100% rejection. As expected, the polarisation
values predicted by our model are only slightly smaller than those obtained using Dresner’s solution. Brian’s finite difference solution [17] for variable wall flux and partial rejection is an extension of Dresner’s Laplace transform solution [19] for constant wall flux and complete rejection. The Leveque solution significantly over-estimates the polarisation because the effects of axial convection are ignored. The integral method includes axial convection but requires knowledge of the shape of the polarisation profile. Taking an exponential profile under-estimates the polarisation while a quadratic profile over-estimates it. The method of rapidly varying boundary conditions (RVBC) with the inclusion of axial convection produces results that are a slight improvement to the latter version of the integral method. The additional effect of variable wall flux (B2 = 0.5 or 0.25) on various solution methods in shown in Fig. 4. The data from our model is close to but below the predictions of the integral method which Doshi et al. [30] conclude to be a satisfactory solution technique for these conditions while the Leveque and RVBC methods become less accurate as B2 increases. At a fixed value of the wall concentration, if B2 is larger, the concentration boundary layer is
132
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Fig. 3. Comparison of the CFD data with the data of Doshi et al. [30] showing the effect of various assumptions on the evaluation of concentration polarisation for constant wall flux, constant solution properties and 100% rejection.
thinner because the wall flux is lower and there is less cumulative extraction of permeate. For thin boundary layers, the shape of the polarisation profile will have less of an effect on the accuracy of results obtained
using the integral method. Thus, our model and the integral method agree well when there are large reductions in the permeation with concentration and, hence, distance, but the agreement diminishes as the
Fig. 4. Comparison of the CFD data with the data of Doshi et al. [30] showing the effect of variable wall flux (B2 = 0.5 or 0.25) on concentration polarisation calculated using different methods for constant solution properties and 100% rejection.
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
133
Fig. 5. Comparison of CFD data with that of Doshi et al. [30] for variable wall flux (B2 = 0.135), variable solution properties (B1 = 0.114) and 100% rejection.
reduction in permeation becomes smaller (i.e. B2 becomes smaller). In Fig. 5, the effect of variable wall flux (B2 = 0.135), variable physical properties (B0 = 0.114) and 100% rejection on concentration polarisation are investigated. To obtain data for this figure, Doshi et al. [30] used the following expression in which viscosity variations are inversely proportional to variations in diffusivity as predicted by the Wilke-Chang correlation: D µ0 1 − B0 (c/c0 ) 2.5 = (8) = D0 µ 1 − B0 It should be noted that, for the example shown here, the results are only valid up to concentrations of approximately 8.5 wt.% where the diffusivity becomes zero and the viscosity becomes infinite. In addition, although reasonable at low concentrations, the Wilke-Change correlation usually loses validity at the higher concentrations (>10 wt.%) typically encountered in industrial membrane processes [35]. For comparison purposes, we initially implemented the physical property variation Eq. (8) in our CFD code. From Fig. 5, it can be seen that our model consistently predicts polarisation values that are much
lower than those of Doshi et al. [30]. As the value of B2 used for the comparison in Fig. 5 is much lower than those used for Fig. 4, the discrepancy between our data and that of Doshi et al. [30] increases as expected for constant properties. When physical property variations are added to the calculations, the discrepancy between the data sets increases remarkably. The RVBC method includes axial variations in viscosity and diffusivity but ignores transverse variations. For thin boundary layers and low values of the wall concentration, the assumption of no transverse variation in physical properties is less important. As wall concentrations increase, changes to the shape of the polarisation profile due to transverse property variations lead to an over-estimation of the polarisation by the RVBC method. However, the effect of the different property variations on concentration polarisation changes in the same order for both our model and the RVBC model. The order is a result of the particular model chosen for the physical properties and may be different when other equations are used to describe the property variations. Further insights into the effect of variations in physical properties on the concentration polarisation can be obtained from Figs. 6 and 7. In these
134
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Fig. 6. Comparison of CFD data with that of Gill et al. [35] for constant wall flux (10−5 m s−1 ), constant or variable solution properties and 100% rejection.
Fig. 7. The effect of variable wall flux (β = 0.25), variable solution properties and 90 or 100% constant rejection values.
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
figures, the variation in viscosity and diffusivity is not confined to the moderate levels defined by the Wilke-Chang correlation but is more representative of the behaviour encountered over the large range of concentrations typically encountered in industrial membrane processes. For the figures, a model solute was defined with the following physical properties for viscosity µ = 0.0009086 e0.00244c
2
(9)
and for diffusivity D=
4.46 × 10−11 tanh(0.159c) c
(10)
Fig. 6 shows the effect of variations in physical properties for a constant wall flux (vw = 10−5 m s−1 ) and 100% rejection. Again, just as for the comparison with the data of Doshi et al. [30], there is a slight discrepancy between the values obtained using our model and those obtained by Gill et al. [35] for constant viscosity. For variable viscosity, the difference between the models increases dramatically. This result is to be expected because the solution developed by Gill et al. [35] uses a variable substitution conversion to translate the equations for variable properties into the same format as that obtained by Dresner [19] for constant properties. Then, Dresner’s constant property solution [19] is integrated numerically to obtain estimates of the polarisation for variable properties. This derivation is acceptable for variable viscosity but produces counter-intuitive values for variable diffusivity because the variable substitution conversion breaks down as the diffusivity variation increases with concentration. The data from our model show trends for the effect of different property variations opposite to that observed for Doshi et al. [30] because of the form of the property variation Eqs. (9) and (10). The additional effects of different values of constant rejection (100 and 90%) and variable wall flux (β = 0.25) on concentration polarisation are shown in Fig. 7. Here it is seen that changes in rejection impose relatively small effects in addition to variable properties on concentration polarisation but variations in flux along the membrane wall significantly affect the polarisation.
135
4. Conclusions In this paper, we have developed a general purpose CFD model of concentration polarisation and fluid flow in pressure driven membrane separation processes. The model has been rigorously tested and validated and shown to perform correctly from a quantitative and qualitative perspective. Now that the model has been validated it can be confidently applied to more complex situations. The validation process has highlighted the need for very fine meshes near the wall and the use of suitable high order numerical schemes, especially when the polarisation is high. In fact, the use of approximate solutions for estimation of concentration polarisation when the physical properties vary with concentration is highly suspect. Correct modelling of rejection and property variations is also essential. This is difficult because there is no a priori way of determining the variations under the wide range of conditions achieved at the membrane wall. While viscosity and diffusivity can be measured in high concentrations solutions, such solutions may not be the same as conditions in the boundary layer which are affected by sieving and/or shear denaturation. While the CFD model developed contains specific equations for the cases investigated, it can be readily modified to encompass any combination of variations in wall flux, rejection, viscosity and diffusivity. The flexibility of the current model is only limited by the ability of the user to accurately define the property variations in real world applications. Acknowledgements One of us (DEW) would like to thank the Department of Chemical Engineering at the University of Sydney, Australia for arranging visiting privileges. References [1] W.F. Blatt, A. Dravid, A.S. Michaels, L. Nelson, Solute polarization and cake formation in membrane ultrafiltration: causes, consequences and control techniques, in: J.E. Flinn (Ed.), Membrane Science and Technology, Plenum Press, New York, 1970, pp. 47–97. [2] M.C. Porter, Concentration polarization with membrane ultrafiltration, Ind. Eng. Chem. Fund. 11 (1972) 234.
136
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
[3] A.S. Michaels, Effects of feed-side polarization on pervaporative stripping of volatile organic solutes from dilutes solutions: a generalized analytical treatment, J. Memb. Sci. 101 (1995) 117. [4] S. Bhattacharya, S.-T. Hwang, Concentration polarization, separation factor, and Peclet number in membrane processes, J. Memb. Sci. 132 (1997) 73. [5] A.L. Zydney, Stagnant film model for concentration polarization in membrane systems, J. Memb. Sci. 130 (1997) 275. [6] S.P. Agashichev, Modelling concentration polarization phenomena for membrane channel with cylindrical geometry in an ultrafiltration process, Desalination 119 (1998) 159. [7] A.S. Berman, Laminar flow in channels with porous walls, J. Appl. Phys. 24 (1953) 1232. [8] S.W. Yuan, A.B. Finkelstein, Laminar flow with injection and suction through a porous wall, Trans. ASME 78 (1956) 719. [9] R.M. Terrill, Laminar flow in a uniformly porous channel, Aeronautical Q. 15 (1964) 299. [10] T.K. Sherwood, P.L.T. Brian, R.E. Fisher, L. Dresner, Salt concentration at phase boundaries in desalination by reverse osmosis, Ind. Eng. Chem. Fund. 4 (1965) 113. [11] M. Friedman, J. Gillis, Viscous flow in a pipe with absorbing walls, J. Appl. Mech. 34 (1967) 819. [12] R.M. Terrill, G. Walker, Heat and mass transfer in laminar flow between parallel porous plates, Appl. Sci. Res. 18 (1967) 193. [13] S. Kimura, S. Sourirajan, Concentration polarization effects in reverse osmosis using porous cellulose acetate membranes, Ind. Eng. Chem. Process Des. Dev. 7 (1968) 41. [14] T. Mizushina, S. Takeshita, G. Unno, Study of flow in a porous tube with radial mass flux, J. Chem. Eng. Japan 4 (1971) 135. [15] L.S. Galowin, L.S.andM.J. DeSantis, Investigation of laminar flow in a porous pipe with variable wall suction, AIAA J. 12 (1974) 1585. [16] J.P. Quaile, E.K. Levy, Laminar flow in a porous tube with suction, J. Heat Transfer 97 (1975) 66. [17] P.L.T. Brian, Concentration polarization in reverse osmosis desalination with variable flux and incomplete salt rejection, Ind. Eng. Chem. Fund. 4 (1965) 439. [18] T. Kotzev, Numerical study of the fluid dynamics and mass transfer of an ultrafiltration performance in a tube membrane module, Int. J. Eng. Sci. 32 (1994) 359. [19] L. Dresner, Boundary Layer Build-Up in the Demineralization of Salt Water by Reverse Osmosis, Oak Ridge Nat. Lab., Report 3621, May 1964. [20] W.N. Gill, C. Tien, D.W. Zeh, Concentration polarization effects in a reverse osmosis system, Ind. Eng. Chem. Fundam. 4 (1965) 433. [21] W.N. Gill, C. Tien, D.W. Zeh, Analysis of continuous reverse osmosis systems for desalination, Int. J. Heat Mass Transfer 9 (1966) 907. [22] A.A. Kozinski, F.P. Schmidt, E.N. Lightfoot, Velocity profiles in porous-walled ducts, Ind. Eng. Chem. Fundam. 9 (1970) 502. [23] R.M. Terrill, Laminar flow in a porous tube, J. Fluids Eng. 105 (1983) 303.
[24] R. Singh, R.L. Laurence, Influence of slip velocity at a membrane surface on ultrafiltration performance. I. Channel flow systems, Int. J. Heat Mass Transfer 22 (1979) 721. [25] H.M. Yeh, T.W. Cheng, Analysis of the slip effect on the permeate flux in membrane ultrafiltration, J. Memb. Sci. 154 (1999) 41. [26] A. Acrivos, On the solution of the convection equation in laminar boundary layer flows, Chem. Eng. Sci. 17 (1961) 457. [27] D.R. Olander, The influence of physical property variations on liquid-phase mass transfer for various laminar flows, Int. J. Heat Mass Transfer 5 (1962) 765. [28] A.S. Emanuel, D.R. Olander, High flux solid–liquid mass transfer, Int. J. Heat Mass Transfer 7 (1964) 539. [29] M.R. Goldman, L.R. Barrett, Mass transfer with concentration-dependent properties in streamline flow, Trans. Inst. Chem. Engrs. 47 (1969) T29. [30] M.R. Doshi, K. Dewan, W.N. Gill, The effect of concentration dependent viscosity and diffusivity on concentration polarization in reverse osmosis flow systems, AIChE Symp. Ser. 68 (1971) 323. [31] A.A. Kozinski, E.N. Lightfoot, Ultrafiltration of proteins in stagnation flow, AIChE J. 17 (1971) 81. [32] Y.T. Shah, Mass transport in reverse osmosis in case of variable diffusivity, Int. J. Heat Mass Transfer 14 (1971) 921. [33] J.J.S. Shen, R.F. Probstein, Ultrafiltration of macromolecular solutions, Ind. Eng. Chem. Fund. 16 (1972) 459. [34] R.G. Wylie, The transfer coefficients of a laminar boundary layer with variable fluid properties, Chem. Eng. J. 6 (1973) 1. [35] W.N. Gill, D.E. Wiley, C.J.D. Fell, A.G. Fane, Effect of viscosity on concentration polarization in ultrafiltration, AIChE J. 34 (1988) 1563. [36] S. De, S. Bhattacharjee, A. Sharma, P.K. Bhattacharya, Generalized integral and similarity solutions of the concentration profiles for osmotic pressure controlled ultrafiltration, J. Memb. Sci. 130 (1997) 99. [37] L. Derzansky, W.N. Gill, The mechanisms of brine-side mass transfer in a horizontal tubular membrane, AIChE J. 20 (1974) 751. [38] K.H. Youm, A.G. Fane, D.E. Wiley, Effects of natural convection instability on membrane performance in dead-end and cross-flow ultrafiltration, J. Memb Sci. 116 (1996) 229. [39] C. Kleinstreuer, G. Belfort, Mathematical modeling of fluid flow and solute distribution in pressure-driven membrane modules, in: G. Belfort (Ed.), Synthetic Membrane Processes—Fundamentals and Water Applications, Academic Press, New York, 1984, pp. 131–190. [40] L. Huang, M.T. Morrissey, Finite element analysis as a tool for crossflow membrane filter simulation, J. Memb. Sci. 155 (1999) 19. [41] Y. Lee, M.M. Clark, A numerical model of steady-state permeate flux during cross-flow ultrafiltration, Desalination 109 (1997) 241. [42] M. Hansen, V.A. Barker, O. Hassager, Spectral element simulation of ultrafiltration, Chem. Eng. Sci. 53 (1998) 3099. [43] V.M. Geraldes, V.A. Semião, M.N. de Pinho, Nanofiltration mass transfer at the entrance region of a slit laminar flow, Ind. Eng. Chem. Res. 37 (1998) 4792.
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137 [44] V. Geraldes, V. Semião, M. Norberta Pinho, Numerical modelling of mass transfer in slits with semi-permeable membrane walls, Eng. Comput. 17 (2000) 192. [45] V. Geraldes, V. Semião, M. Norberta de Pinho, Flow and mass transfer modelling in nanofiltration, J. Memb. Sci. 191 (2001) 109. [46] S.K. Karode, A. Kumar, Flow visualization through spacer filled channels by computational fluid dynamics. I. Pressure
137
drop and shear rate calculations for flat sheet geometry, J. Memb. Sci. 193 (2001) 69. [47] J.P. Van Doormal, G.D. Raithby, Enhancements of the SIMPLE method for predicting incompressible fluid flows, Numer. Heat transfer 7 (1984) 147. [48] C. Hirsch, Numerical Computation of Internal and External Flows, vol. 2, Wiley, New York, 1990, pp. 536–545.
Techniques for computational fluid dynamics modelling of flow in membrane channels Dianne E. Wiley a,∗ , David F. Fletcher b a
UNESCO Centre for Membrane Science and Technology, School of Chemical Engineering and Industrial Chemistry, University of New South Wales, Sydney 2052, NSW, Australia b Department of Chemical Engineering, University of Sydney, Sydney 2006, NSW, Australia Received 30 December 2001; received in revised form 15 May 2002; accepted 5 September 2002
Abstract Accurate modelling of the flow and concentration polarisation in pressure driven membrane processes is inhibited by the complex couplings in the flow equations along with any added effects of variable solution properties. A generic computational fluid dynamics (CFD) model has been developed which incorporates these effects and describes the flow across the membrane wall. The results have been validated against classical solutions available in the literature. Extended work indicates that overly simplified expressions for the dependence of viscosity and diffusivity on concentration produce velocity and concentration profiles that may grossly misrepresent reality. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Concentration polarisation; CFD; Theory; Ultrafiltration; Reverse osmosis
1. Introduction The modelling of flow and concentration polarisation in membrane or thin channels with permeable walls is not new. Numerous approximate one-dimensional [1–6] and two-dimensional [7–38] models of varying degrees of complexity and simplification have appeared in the literature. A succinct review of much of the early work was prepared by Kleinstreuer and Belfort [39]. The attraction of these analytical and semi-analytical models is the ability to tailor them for relatively quick and easy investigations of the effect of different approximations and parameters on the performance of membrane ∗ Corresponding author. Tel.: +61-2-98385-4304; fax: +61-2-98385-5966. E-mail address: [email protected] (D.E. Wiley).
systems. Models exist which permit investigation of the effects of variable permeation [17–23], changes in rejection [17,18], effects of wall slip [24,25], variable physical properties of the solution being processed [26–36] and gravitational effects [37,38]. Despite the simplifications, many of the models have proved to be remarkably accurate at predicting membrane performance under some limited circumstances. Development of a generic model with widespread applicability to membrane systems requires much more rigorous and robust techniques than those described above. Computational fluid dynamics (CFD) provides one such technique. There have been many fewer attempts at modelling membrane systems in this manner. This is most likely due to the complexity of the problem and the fact that the ODE based approaches described previously are computationally very efficient. However, CFD models are needed if
0376-7388/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 6 - 7 3 8 8 ( 0 2 ) 0 0 4 1 2 - X
128
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Nomenclature B0 B1 B2 c C+ D h H L p R Re Sc u v x X y
constant in Eq. (8) constant in Eq. (5) constant in Eq. (5) concentration (wt.%) dimensionless concentration, c/c0 diffusivity (m2 s−1 ) half flow channel height (m) half permeate channel height (m) channel length (m) pressure (Pa) rejection coefficient Reynolds number Schmidt number velocity in the axial direction (m s−1 ) velocity in the transverse direction (m s−1 ) axial direction (m) fractional water removal, vw0 x/u0 h transverse direction (m)
Greek letters α D/vw h β constant in Eq. (4) µ dynamic viscosity (Pa s) ν kinematic viscosity (m2 s−1 ) ρ density (kg m−3 ) ξ 3B13 x/ h Sc Re Subscripts p permeate w wall 0 inlet
the flow in complex geometries, which include recirculation, are to be studied. A variety of techniques have been used and are summarised below. Huang and Morrissey [40] used the finite element technique to solve for the concentration field given an approximate solution to the velocity field, by neglecting stream-wise diffusion and assuming constant properties. This work is very limited but it highlighted the need for extremely high resolution near the membrane surface because of the extremely thin concentration boundary layer. Lee and Clark [41] performed a similar analysis but this time solved the concentration equation using a finite difference method. They also
included a shear induced diffusion term, as well as the usual molecular diffusion term, in the concentration equation. In addition, they developed an expression for the permeate flux which is based on diffusion through a stagnant cake layer assumed to develop on the membrane surface. Both works are limited by the fact that they assume the velocity field to be known. Hansen et al. [42] used a spectral method to solve the laminar Navier-Stokes equations for channel flow with permeate flow through a wall. They allowed the permeate velocity to be a function of the static and osmotic pressure differences and made the viscosity depend on concentration. The spectral method used had high spatial accuracy and required relatively few elements but used polynomials of degree five for concentration and pressure and three for velocity, giving a large number of degrees of freedom. Whilst the method worked well for the parabolic flow considered and for the case of 100% rejection, it is not clear how easily it can be extended to complex geometries and to cases where the rejection is less than 100%. Geraldes et al. [43,44] performed finite volume simulations for permeation of water through a membrane for a simple channel flow. They solved both the flow and the concentration fields for a constant property situation. The simulations used a solute of salt, so that the diffusivity was relatively large. Despite this they needed very fine meshes to obtain a grid independent solution. They present considerable discussion on differencing schemes but their conclusions are of limited value, as they apply to a flow-aligned geometry with a relatively thick concentration boundary layer. Later work [45] has extended this model to include variable fluid properties. The model has been validated against experimental data on nanofiltration and used to fit effective rejection coefficients. There are almost no simulations of flow in thin channels or membrane systems in three dimensions. Recent work by Karode and Kumar [46] using relatively coarse grids and first order differencing for convection has investigated the effect of spacers on pressure drop and shear stress, but as yet does not include the effect of permeation on flow and concentration polarisation. The present work is directed towards removing some of the assumptions discussed above and to validate a general purpose finite volume model for use in the simulation of membrane systems. As a first step, we performed simulations in two dimensions
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
to enable validation of the modelling framework but note that this approach is easily extended to certain three-dimensional flows.
2. Theory The starting point for the work presented here is the commercially-available finite volume code CFX4 (CFX, 1998) which solves the conservation equations on a co-located grid using the SIMPLEC algorithm [47]. Specifically, we have assumed the flow to be steady, incompressible with constant density, and laminar. The transport properties µ (viscosity) and D (diffusivity) can be a function of the concentration. Given these assumptions, the equations for conservation of mass, momentum and solute concentration are: ∂u ∂v + =0 ∂x ∂y
(1)
∂u ∂u ρ u +v ∂x ∂y ∂p ∂ ∂u ∂ ∂u =− + µ + µ ∂x ∂x ∂x ∂y ∂y ∂v ∂v ρ u +v ∂x ∂y ∂p ∂ ∂v ∂ ∂v =− + µ + µ ∂y ∂x ∂x ∂y ∂y
When the concentration polarisation is large, it is important that this condition is applied well downstream of the end of the membrane, otherwise the solution over the last part of the membrane can be affected significantly. At the symmetry plane, v = 0, and the normal gradients of the tangential velocity (u) and the concentration (c) are set to zero. At the walls, the no slip condition is imposed and the normal derivative of the concentration is set to zero. At a membrane, the conditions are more complex, as flow permeates through the wall. For a membrane located at the plane y = 0, the tangential velocity (u) is set to zero and the normal velocity is set to vw . In some cases the wall (permeate) velocity was set to a constant value and in others it was made a function of the wall concentration of solute. The expression used was either that of Brian [17], which takes the form cw vw = v0 1 − β (4) c0 or that of Doshi et al. [30] in which v0 B1 cw 1 − B2 vw = h Sc c0
(2)
and ∂c ∂ ∂c ∂ ∂c ∂c = D + D (3) u +v ∂x ∂y ∂x ∂x ∂y ∂y 2.1. Boundary conditions
129
(5)
The boundary condition on the concentration results from a balance of the convective and diffusive fluxes. In addition, it must take into account the fact that not all of the solute permeates through the membrane. This is done via the use of a rejection coefficient, R, and the concentration boundary condition is given by D
∂c + vw Rcw = 0 ∂y
on the feed side of the membrane and ∂c D + vw cp = (1 − R)cw vw ∂y
(6)
(7)
on the permeate side of the membrane. At the inlet, the flow is assumed to be fully developed and a parabolic flow is specified. As this boundary condition is applied at the start of the membrane, then a transverse velocity, as well as a normal velocity is specified, based on the solution of Berman [7]. Other arbitrary inlet flow conditions could be specified for different applications with appropriate modification of the inlet channel length if required. A uniform inflow concentration of c0 is specified. At the outlet, all derivatives in the flow direction are set to zero.
2.2. The geometry and computational mesh used A model, shown schematically in Fig. 1, was set up to represent a section of a channel with a membrane separating two regions. The physical dimensions of the geometry were varied to match the conditions in the cited references. For example, in the case of the comparison with the data of Brian [17] the channel was 2.7 m long, the feed channel had a half channel
130
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Fig. 1. Geometrical model of the flow channel used for the simulations.
height of 1 mm and the permeate channel had a height of 0.6 mm. Flow is imposed in the top section of the channel and occurs in the bottom channel as mass is transferred through the membrane. Typically a mesh with 50 cells across the half channel, 10 cells across the permeate channel and 60 cells along the channel were used. A non-uniform mesh was used, with the mesh density being higher at the start of the membrane and near to the surface of the membrane. Care was taken to avoid excessively high cell stretch rates. In all cases, we checked the grid independence of the solution and made sure that there were sufficient cells in the concentration polarisation layer that the thickness of this layer and the peak concentration were independent of the mesh used. 2.3. Numerical procedure All convective terms were differenced using a MUSCL scheme with the van Leer limiter to provide high spatial accuracy without the introduction of unphysical dispersive errors [48]. Standard second order differencing was used for the diffusive terms. Stone’s method was used to solve the resulting matrix equations for the velocity and concentration fields and incomplete Cholesky conjugate gradient (ICCG) was used for the pressure. Under-relaxation factors of 0.65 on the velocity components and unity on the pressure and concentration worked well for situations in which the polarisation was not great, such as those simulated by Brian [17]. However, when very high concentrations developed near the wall, resulting in a very large fluid viscosity, it was necessary to reduce the relaxation factors on all but the pressure equation to account for the extreme non-linear coupling
between the viscosity and the solute concentration. In extreme cases, values as low as 0.05 were required, resulting in the need to perform several thousand iterations to obtain convergence. In some circumstances we also found it to be beneficial to under-relax the viscosity. However, in all cases we were able to obtain extremely good convergence, with the total mass residual being reduced to 0.05% of the inflow mass.
3. Results and discussion The model developed in this paper has been validated against a number of ‘classical’ semi-analytical solutions that are well accepted in the membrane literature. Comparisons to be shown encompass the effects of various combinations of variable and constant wall flux, variable and constant solution properties and several values of constant rejection. Unfortunately there is no single solution available in the literature that covers all of the combinations. A range of validations of increasing complexity has therefore been selected to demonstrate the general validity and utility of the model developed here. In comparing the data, it must be remembered that the CFD model performs a rigorous evaluation of the complete transport equations presented in Eqs. (1)–(3) while the semi-analytical methods each approximate the solution of the equations in different ways. Thus, some discrepancies between the various methods would be anticipated. This will be elaborated on in the ensuing discussion. In order to demonstrate that the boundary conditions at the membrane have been correctly specified, Fig. 2 compares our CFD solutions with data published by Brian [17] for variable wall flux (α = 0.27, β = 1),
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
131
Fig. 2. Comparison of the CFD data with the data of Brian for variable wall flux (α = 0.27, β = 1), constant solution properties and constant rejection (100 or 90%).
constant physical properties of the solution (ρ, D, µ) and constant rejections of 100 or 90%. Brian [17] used a finite difference technique with truncated versions of the velocity profiles determined by Berman [7] using a perturbation technique. Both Brian [17] and Berman [7] demonstrated the acceptability of truncating the velocity profiles when only small amounts of fluid were removed over the length of the permeating channel. Fig. 2 shows the excellent numerical accuracy of the model developed, although, it consistently predicts slightly smaller values of the wall concentration. The CFD model solves the full Navier-Stokes and concentration equations, whereas Brian [17] uses a truncated velocity distribution and ignores stream-wise diffusion of the solute. Comparison of the calculated velocity profile with the truncated expression shows a difference of less than three percent in the wall region for the cases investigated in this paper. This difference, together with the neglect of stream-wise diffusion, appears sufficient to explain the minor differences. Fig. 3 compares the effect of some approximate analytical solutions presented by Doshi et al. [30] against the rigorous CFD solution for constant wall flux (B2 = 0), constant physical properties (B0 = 0) and 100% rejection. As expected, the polarisation
values predicted by our model are only slightly smaller than those obtained using Dresner’s solution. Brian’s finite difference solution [17] for variable wall flux and partial rejection is an extension of Dresner’s Laplace transform solution [19] for constant wall flux and complete rejection. The Leveque solution significantly over-estimates the polarisation because the effects of axial convection are ignored. The integral method includes axial convection but requires knowledge of the shape of the polarisation profile. Taking an exponential profile under-estimates the polarisation while a quadratic profile over-estimates it. The method of rapidly varying boundary conditions (RVBC) with the inclusion of axial convection produces results that are a slight improvement to the latter version of the integral method. The additional effect of variable wall flux (B2 = 0.5 or 0.25) on various solution methods in shown in Fig. 4. The data from our model is close to but below the predictions of the integral method which Doshi et al. [30] conclude to be a satisfactory solution technique for these conditions while the Leveque and RVBC methods become less accurate as B2 increases. At a fixed value of the wall concentration, if B2 is larger, the concentration boundary layer is
132
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Fig. 3. Comparison of the CFD data with the data of Doshi et al. [30] showing the effect of various assumptions on the evaluation of concentration polarisation for constant wall flux, constant solution properties and 100% rejection.
thinner because the wall flux is lower and there is less cumulative extraction of permeate. For thin boundary layers, the shape of the polarisation profile will have less of an effect on the accuracy of results obtained
using the integral method. Thus, our model and the integral method agree well when there are large reductions in the permeation with concentration and, hence, distance, but the agreement diminishes as the
Fig. 4. Comparison of the CFD data with the data of Doshi et al. [30] showing the effect of variable wall flux (B2 = 0.5 or 0.25) on concentration polarisation calculated using different methods for constant solution properties and 100% rejection.
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
133
Fig. 5. Comparison of CFD data with that of Doshi et al. [30] for variable wall flux (B2 = 0.135), variable solution properties (B1 = 0.114) and 100% rejection.
reduction in permeation becomes smaller (i.e. B2 becomes smaller). In Fig. 5, the effect of variable wall flux (B2 = 0.135), variable physical properties (B0 = 0.114) and 100% rejection on concentration polarisation are investigated. To obtain data for this figure, Doshi et al. [30] used the following expression in which viscosity variations are inversely proportional to variations in diffusivity as predicted by the Wilke-Chang correlation: D µ0 1 − B0 (c/c0 ) 2.5 = (8) = D0 µ 1 − B0 It should be noted that, for the example shown here, the results are only valid up to concentrations of approximately 8.5 wt.% where the diffusivity becomes zero and the viscosity becomes infinite. In addition, although reasonable at low concentrations, the Wilke-Change correlation usually loses validity at the higher concentrations (>10 wt.%) typically encountered in industrial membrane processes [35]. For comparison purposes, we initially implemented the physical property variation Eq. (8) in our CFD code. From Fig. 5, it can be seen that our model consistently predicts polarisation values that are much
lower than those of Doshi et al. [30]. As the value of B2 used for the comparison in Fig. 5 is much lower than those used for Fig. 4, the discrepancy between our data and that of Doshi et al. [30] increases as expected for constant properties. When physical property variations are added to the calculations, the discrepancy between the data sets increases remarkably. The RVBC method includes axial variations in viscosity and diffusivity but ignores transverse variations. For thin boundary layers and low values of the wall concentration, the assumption of no transverse variation in physical properties is less important. As wall concentrations increase, changes to the shape of the polarisation profile due to transverse property variations lead to an over-estimation of the polarisation by the RVBC method. However, the effect of the different property variations on concentration polarisation changes in the same order for both our model and the RVBC model. The order is a result of the particular model chosen for the physical properties and may be different when other equations are used to describe the property variations. Further insights into the effect of variations in physical properties on the concentration polarisation can be obtained from Figs. 6 and 7. In these
134
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
Fig. 6. Comparison of CFD data with that of Gill et al. [35] for constant wall flux (10−5 m s−1 ), constant or variable solution properties and 100% rejection.
Fig. 7. The effect of variable wall flux (β = 0.25), variable solution properties and 90 or 100% constant rejection values.
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
figures, the variation in viscosity and diffusivity is not confined to the moderate levels defined by the Wilke-Chang correlation but is more representative of the behaviour encountered over the large range of concentrations typically encountered in industrial membrane processes. For the figures, a model solute was defined with the following physical properties for viscosity µ = 0.0009086 e0.00244c
2
(9)
and for diffusivity D=
4.46 × 10−11 tanh(0.159c) c
(10)
Fig. 6 shows the effect of variations in physical properties for a constant wall flux (vw = 10−5 m s−1 ) and 100% rejection. Again, just as for the comparison with the data of Doshi et al. [30], there is a slight discrepancy between the values obtained using our model and those obtained by Gill et al. [35] for constant viscosity. For variable viscosity, the difference between the models increases dramatically. This result is to be expected because the solution developed by Gill et al. [35] uses a variable substitution conversion to translate the equations for variable properties into the same format as that obtained by Dresner [19] for constant properties. Then, Dresner’s constant property solution [19] is integrated numerically to obtain estimates of the polarisation for variable properties. This derivation is acceptable for variable viscosity but produces counter-intuitive values for variable diffusivity because the variable substitution conversion breaks down as the diffusivity variation increases with concentration. The data from our model show trends for the effect of different property variations opposite to that observed for Doshi et al. [30] because of the form of the property variation Eqs. (9) and (10). The additional effects of different values of constant rejection (100 and 90%) and variable wall flux (β = 0.25) on concentration polarisation are shown in Fig. 7. Here it is seen that changes in rejection impose relatively small effects in addition to variable properties on concentration polarisation but variations in flux along the membrane wall significantly affect the polarisation.
135
4. Conclusions In this paper, we have developed a general purpose CFD model of concentration polarisation and fluid flow in pressure driven membrane separation processes. The model has been rigorously tested and validated and shown to perform correctly from a quantitative and qualitative perspective. Now that the model has been validated it can be confidently applied to more complex situations. The validation process has highlighted the need for very fine meshes near the wall and the use of suitable high order numerical schemes, especially when the polarisation is high. In fact, the use of approximate solutions for estimation of concentration polarisation when the physical properties vary with concentration is highly suspect. Correct modelling of rejection and property variations is also essential. This is difficult because there is no a priori way of determining the variations under the wide range of conditions achieved at the membrane wall. While viscosity and diffusivity can be measured in high concentrations solutions, such solutions may not be the same as conditions in the boundary layer which are affected by sieving and/or shear denaturation. While the CFD model developed contains specific equations for the cases investigated, it can be readily modified to encompass any combination of variations in wall flux, rejection, viscosity and diffusivity. The flexibility of the current model is only limited by the ability of the user to accurately define the property variations in real world applications. Acknowledgements One of us (DEW) would like to thank the Department of Chemical Engineering at the University of Sydney, Australia for arranging visiting privileges. References [1] W.F. Blatt, A. Dravid, A.S. Michaels, L. Nelson, Solute polarization and cake formation in membrane ultrafiltration: causes, consequences and control techniques, in: J.E. Flinn (Ed.), Membrane Science and Technology, Plenum Press, New York, 1970, pp. 47–97. [2] M.C. Porter, Concentration polarization with membrane ultrafiltration, Ind. Eng. Chem. Fund. 11 (1972) 234.
136
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137
[3] A.S. Michaels, Effects of feed-side polarization on pervaporative stripping of volatile organic solutes from dilutes solutions: a generalized analytical treatment, J. Memb. Sci. 101 (1995) 117. [4] S. Bhattacharya, S.-T. Hwang, Concentration polarization, separation factor, and Peclet number in membrane processes, J. Memb. Sci. 132 (1997) 73. [5] A.L. Zydney, Stagnant film model for concentration polarization in membrane systems, J. Memb. Sci. 130 (1997) 275. [6] S.P. Agashichev, Modelling concentration polarization phenomena for membrane channel with cylindrical geometry in an ultrafiltration process, Desalination 119 (1998) 159. [7] A.S. Berman, Laminar flow in channels with porous walls, J. Appl. Phys. 24 (1953) 1232. [8] S.W. Yuan, A.B. Finkelstein, Laminar flow with injection and suction through a porous wall, Trans. ASME 78 (1956) 719. [9] R.M. Terrill, Laminar flow in a uniformly porous channel, Aeronautical Q. 15 (1964) 299. [10] T.K. Sherwood, P.L.T. Brian, R.E. Fisher, L. Dresner, Salt concentration at phase boundaries in desalination by reverse osmosis, Ind. Eng. Chem. Fund. 4 (1965) 113. [11] M. Friedman, J. Gillis, Viscous flow in a pipe with absorbing walls, J. Appl. Mech. 34 (1967) 819. [12] R.M. Terrill, G. Walker, Heat and mass transfer in laminar flow between parallel porous plates, Appl. Sci. Res. 18 (1967) 193. [13] S. Kimura, S. Sourirajan, Concentration polarization effects in reverse osmosis using porous cellulose acetate membranes, Ind. Eng. Chem. Process Des. Dev. 7 (1968) 41. [14] T. Mizushina, S. Takeshita, G. Unno, Study of flow in a porous tube with radial mass flux, J. Chem. Eng. Japan 4 (1971) 135. [15] L.S. Galowin, L.S.andM.J. DeSantis, Investigation of laminar flow in a porous pipe with variable wall suction, AIAA J. 12 (1974) 1585. [16] J.P. Quaile, E.K. Levy, Laminar flow in a porous tube with suction, J. Heat Transfer 97 (1975) 66. [17] P.L.T. Brian, Concentration polarization in reverse osmosis desalination with variable flux and incomplete salt rejection, Ind. Eng. Chem. Fund. 4 (1965) 439. [18] T. Kotzev, Numerical study of the fluid dynamics and mass transfer of an ultrafiltration performance in a tube membrane module, Int. J. Eng. Sci. 32 (1994) 359. [19] L. Dresner, Boundary Layer Build-Up in the Demineralization of Salt Water by Reverse Osmosis, Oak Ridge Nat. Lab., Report 3621, May 1964. [20] W.N. Gill, C. Tien, D.W. Zeh, Concentration polarization effects in a reverse osmosis system, Ind. Eng. Chem. Fundam. 4 (1965) 433. [21] W.N. Gill, C. Tien, D.W. Zeh, Analysis of continuous reverse osmosis systems for desalination, Int. J. Heat Mass Transfer 9 (1966) 907. [22] A.A. Kozinski, F.P. Schmidt, E.N. Lightfoot, Velocity profiles in porous-walled ducts, Ind. Eng. Chem. Fundam. 9 (1970) 502. [23] R.M. Terrill, Laminar flow in a porous tube, J. Fluids Eng. 105 (1983) 303.
[24] R. Singh, R.L. Laurence, Influence of slip velocity at a membrane surface on ultrafiltration performance. I. Channel flow systems, Int. J. Heat Mass Transfer 22 (1979) 721. [25] H.M. Yeh, T.W. Cheng, Analysis of the slip effect on the permeate flux in membrane ultrafiltration, J. Memb. Sci. 154 (1999) 41. [26] A. Acrivos, On the solution of the convection equation in laminar boundary layer flows, Chem. Eng. Sci. 17 (1961) 457. [27] D.R. Olander, The influence of physical property variations on liquid-phase mass transfer for various laminar flows, Int. J. Heat Mass Transfer 5 (1962) 765. [28] A.S. Emanuel, D.R. Olander, High flux solid–liquid mass transfer, Int. J. Heat Mass Transfer 7 (1964) 539. [29] M.R. Goldman, L.R. Barrett, Mass transfer with concentration-dependent properties in streamline flow, Trans. Inst. Chem. Engrs. 47 (1969) T29. [30] M.R. Doshi, K. Dewan, W.N. Gill, The effect of concentration dependent viscosity and diffusivity on concentration polarization in reverse osmosis flow systems, AIChE Symp. Ser. 68 (1971) 323. [31] A.A. Kozinski, E.N. Lightfoot, Ultrafiltration of proteins in stagnation flow, AIChE J. 17 (1971) 81. [32] Y.T. Shah, Mass transport in reverse osmosis in case of variable diffusivity, Int. J. Heat Mass Transfer 14 (1971) 921. [33] J.J.S. Shen, R.F. Probstein, Ultrafiltration of macromolecular solutions, Ind. Eng. Chem. Fund. 16 (1972) 459. [34] R.G. Wylie, The transfer coefficients of a laminar boundary layer with variable fluid properties, Chem. Eng. J. 6 (1973) 1. [35] W.N. Gill, D.E. Wiley, C.J.D. Fell, A.G. Fane, Effect of viscosity on concentration polarization in ultrafiltration, AIChE J. 34 (1988) 1563. [36] S. De, S. Bhattacharjee, A. Sharma, P.K. Bhattacharya, Generalized integral and similarity solutions of the concentration profiles for osmotic pressure controlled ultrafiltration, J. Memb. Sci. 130 (1997) 99. [37] L. Derzansky, W.N. Gill, The mechanisms of brine-side mass transfer in a horizontal tubular membrane, AIChE J. 20 (1974) 751. [38] K.H. Youm, A.G. Fane, D.E. Wiley, Effects of natural convection instability on membrane performance in dead-end and cross-flow ultrafiltration, J. Memb Sci. 116 (1996) 229. [39] C. Kleinstreuer, G. Belfort, Mathematical modeling of fluid flow and solute distribution in pressure-driven membrane modules, in: G. Belfort (Ed.), Synthetic Membrane Processes—Fundamentals and Water Applications, Academic Press, New York, 1984, pp. 131–190. [40] L. Huang, M.T. Morrissey, Finite element analysis as a tool for crossflow membrane filter simulation, J. Memb. Sci. 155 (1999) 19. [41] Y. Lee, M.M. Clark, A numerical model of steady-state permeate flux during cross-flow ultrafiltration, Desalination 109 (1997) 241. [42] M. Hansen, V.A. Barker, O. Hassager, Spectral element simulation of ultrafiltration, Chem. Eng. Sci. 53 (1998) 3099. [43] V.M. Geraldes, V.A. Semião, M.N. de Pinho, Nanofiltration mass transfer at the entrance region of a slit laminar flow, Ind. Eng. Chem. Res. 37 (1998) 4792.
D.E. Wiley, D.F. Fletcher / Journal of Membrane Science 211 (2003) 127–137 [44] V. Geraldes, V. Semião, M. Norberta Pinho, Numerical modelling of mass transfer in slits with semi-permeable membrane walls, Eng. Comput. 17 (2000) 192. [45] V. Geraldes, V. Semião, M. Norberta de Pinho, Flow and mass transfer modelling in nanofiltration, J. Memb. Sci. 191 (2001) 109. [46] S.K. Karode, A. Kumar, Flow visualization through spacer filled channels by computational fluid dynamics. I. Pressure
137
drop and shear rate calculations for flat sheet geometry, J. Memb. Sci. 193 (2001) 69. [47] J.P. Van Doormal, G.D. Raithby, Enhancements of the SIMPLE method for predicting incompressible fluid flows, Numer. Heat transfer 7 (1984) 147. [48] C. Hirsch, Numerical Computation of Internal and External Flows, vol. 2, Wiley, New York, 1990, pp. 536–545.