Next Article in Journal
Effect of Nanopatterning on Concentration Polarization during Nanofiltration
Next Article in Special Issue
Numerical Modeling in Membrane Processes
Previous Article in Journal
Solid-Contact Potentiometric Anion Sensing Based on Classic Silver/Silver Insoluble Salts Electrodes without Ion-Selective Membrane
Previous Article in Special Issue
A Novel Numerical Procedure to Estimate the Electric Charge in the Pore from Filtration of Single-Salt Solutions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Geometrical Influence on Particle Transport in Cross-Flow Ultrafiltration: Cylindrical and Flat Sheet Membranes

Institute of Biological Information Processing (IBI-4), Forschungszentrum Juelich GmbH, 52425 Jülich, Germany
*
Author to whom correspondence should be addressed.
Membranes 2021, 11(12), 960; https://doi.org/10.3390/membranes11120960
Submission received: 8 November 2021 / Revised: 29 November 2021 / Accepted: 1 December 2021 / Published: 6 December 2021
(This article belongs to the Special Issue Numerical Modelling in Membrane Processes)

Abstract

:
Cross-flow membrane ultrafiltration (UF) is used for the enrichment and purification of small colloidal particles and proteins. We explore the influence of different membrane geometries on the particle transport in, and the efficiency of, inside-out cross-flow UF. For this purpose, we generalize the accurate and numerically efficient modified boundary layer approximation (mBLA) method, developed in recent work by us for a hollow cylindrical membrane, to parallel flat sheet geometries with one or two solvent-permeable membrane sheets. Considering a reference dispersion of Brownian hard spheres where accurate expressions for its transport properties are available, the generalized mBLA method is used to analyze how particle transport and global UF process indicators are affected by varying operating parameters and the membrane geometry. We show that global process indicators including the mean permeate flux, the solvent recovery indicator, and the concentration factor are strongly dependent on the membrane geometry. A key finding is that irrespective of the many input parameters characterizing an UF experiment and its membrane geometry, the process indicators are determined by three independent dimensionless variables only. This finding can be very useful in the design, optimization, and scale-up of UF processes.

1. Introduction

Membrane ultrafiltration (UF) is used in daily life, e.g., in water purification, blood treatment by (artificial) kidneys, and protein enrichment [1,2,3]. It is the pressure-driven membrane filtration of smaller, sub-micron sized particles dispersed in a low-molecular solvent (mostly water) under conditions of larger trans-membrane pressure (TMP) values. Owing to strong Brownian motion of the particles, the dispersion remains practically in thermodynamic equilibrium during UF, and there is a significant osmotic pressure buildup along the membrane surface. This distinguishes UF from microfiltration of larger, typically micron-sized particles where Brownian motion and (equilibrium) osmotic pressure effects are negligible. Instead, non-equilibrium micro-hydrodynamic effects including shear-induced collective diffusion are of importance in microfiltration.
UF is often performed using a cross-flow setup, where a feed dispersion is steadily pumped through a parallel bundle of membrane modules having inlet and outlet ports. The particles-enriched dispersion is collected at the outlet ports. Driven by the TMP, there is a non-homogeneous, particles-enriched diffuse layer formed near the surface of a membrane retaining the particles. This concentration-polarization (CP) layer becomes, in general, more pronounced with increasing distance from the inlet port. The CP layer is determined by the balance of gradient (collective) diffusion of particles away and by the flow advection towards the membrane surface. The enlarged dispersion viscosity inside the CP layer further enhances polarization. The osmotic pressure buildup caused by the CP layer counteracts the applied TMP, which in turn lowers the filtration efficiency.
In addition to CP layer formation, there are unwarranted membrane fouling effects, which lower the filtration efficiency. Fouling is caused by specific physico-chemical interactions between particles and membranes, implying, e.g., stagnant cake layer formation, particle adsorption at the membrane wall, or clogging of the membrane pores [3,4,5,6,7]. The present work focuses on the influence of different membrane geometries on particle transport and generic CP layer formation in UF. We are not dealing here with specific fouling effects. The considered UF operating conditions are such that reversible cake layer formation due to crystallization or jamming is avoided.
Among the different membrane geometries encountered in UF setups, there are two common ones depicted in Figure 1, which are discussed in the present work. The first one is the standard cylindrical shape of a hollow fiber or tubular membrane having circular cross-section, characterized by an inner membrane radius R and axial length L with R L . This is referred to in what follows as a cylindrical membrane (CM). The second one is a flat sheets geometry with a rectangular cross-section of height 2 R and width W, where R W and R L . Here, L denotes again the axial length of the membrane. Regarding the flat sheets geometry, we differentiate between two types of membranes. The first type consists of two parallel flat membrane sheets of area W · L each, referred to as the FMM membrane. The second one has a top membrane sheet and a non-permeable bottom sheet referred to as the substrate sheet (c.f. Figure 1). We refer to this asymmetric flat sheets membrane system as an FMS membrane system.
The CM membrane is distinguished from the FMM membrane by a larger cross-sectional circumference-to-area ratio of 2 / R , while this ratio was equal to 1 / R in the FMM case. For otherwise equal membrane properties, i.e., equal hydraulic permeability L p , R, L, and for equal TMP and mean inlet velocity, this difference in the ratio implies a larger solvent recovery for the CM membrane. In spite of being less fouling-resistant than CM membranes, FMM and FMS membrane systems are often used as plate-and-frame setups in industrial applications, owing to their compact design and being less expensive. Moreover, flat sheet membranes can be operated at larger TMP, since they are commonly supported mechanically by plates or rigid porous spacers on the permeate side [8,9,10]. The FMS module with a transparent substrate sheet is also useful for monitoring flow behavior and fouling [11,12,13,14].
In this work, we analyze differences and similarities in particle transport and flow efficiency of CM, FMM, and FMS membrane systems under UF conditions. For simplicity, we assume that the membranes are ideally particle retentive and that the particles are mono-disperse colloidal hard spheres. For fluid-like hard-sphere dispersions, accurate analytic expressions are available for the equilibrium gradient-diffusion coefficient and viscosity, and for the osmotic pressure. To calculate steady-state flow properties and particle concentration profiles, we use a versatile, semi-analytic modified boundary-layer approximation (mBLA) method, developed in earlier work by us for CM membranes [15], and we generalize it to FMM and FMS geometries. It was shown in [15] that this method provides concentration and flow profiles in excellent agreement with the according results obtained from elaborate finite-element calculations. In the mBLA method, the flow outside the CP boundary layer is determined basically by the, in general, hyperbolic axial pressure profile caused by the permeate flux through the membrane, as described in [16,17,18]. The inner flow solution is obtained in a similar way as in classical film theory [4,5,19,20], however with the concentration-dependence of the dispersion properties accounted for. The mBLA method is computationally fast, and different from computational fluid dynamics methods, it offers analytic insight into the functional behavior of concentration and flow properties. In this context, we note that approximate solutions for pure solvent flow in a FMM channel are discussed in [21,22] and for dispersion flow in [23,24,25].
For varied operating conditions, we investigate how the CM, FMM, and FMS geometries affect the concentration, dispersion flux profiles, and global process indicators, including the mean permeate flux and the solvent recovery indicator. The mBLA results recover, in particular, the experimentally observed dependence of the mean permeate flux on the logarithm of the feed concentrations [9]. A major finding of this study is the dependence of the solvent-recovery indicator on a minimal set of only three input variables, which applies to all considered membrane geometries.
The article is organized as follows. In Section 2, we summarize analytic expressions for the suspension transport properties and osmotic pressure, used as input to the generalized mBLA method. The employed cross-flow UF model and its underlying transport equations are discussed in Section 3. Furthermore, the respective boundary conditions for the CM, FMM, and FMS geometries are explained. Section 4 gives the essentials of the mBLA method generalized to the FMM and FMS geometries, and it includes a discussion of its pros and cons. Our results are presented in Section 5. Results for the CP layer profile and particle transport for fixed TMP are discussed in Section 5.1. Global-process indicators are analyzed in Section 5.2 for varying TMP. The minimal set of the three variables fully characterizing the global process indicators is presented in Section 5.3 and discussed in conjunction with Appendix A. Our conclusions are contained in Section 6.

2. Dispersion of Brownian Hard Spheres

The scope of this study is to analyze the influence of three different membrane geometries on generic CP layer effects and the UF efficiency. Since we were not concerned here with effects on UF arising from specific particle interactions, as a reference dispersion, we employ the model of mono-disperse Brownian hard spheres whose thermodynamic and transport properties are characterized by two parameters only, namely, the volume fraction ϕ = ( 4 π / 3 ) n a 3 and the hard-core radius a. Here, n is the number density of spheres.
The flow conditions considered in this work are such that the single-particle Pèclet number, P e a , obeys [26,27]
P e a γ ˙ * a 2 D 0 1 ,
where γ ˙ * is a characteristic shear rate, taken as the shear rate of the velocity field at the inner membrane wall right at the inlet. Furthermore, D 0 = k B T / ( 6 π η s a ) is the single-particle Stokes–Einstein diffusion coefficient for stick hydrodynamic boundary conditions [26], and η s is the viscosity of the Newtonian solvent. For small P e a , the dispersion is only slightly perturbed away from equilibrium (see Figure 2a), with the Brownian forces acting on the particles dominating the viscous ones. This characterizes the UF region, marked by the shaded area in Figure 2a. To model UF, one is allowed to use the equilibrium forms of the transport coefficients and osmotic pressure for a quiescent dispersion, for which accurate analytic expressions are available.
The equilibrium osmotic pressure, Π , of Brownian hard spheres can be described by the Carnahan–Starling expression [28,29]
Π ( ϕ , a ) = k B T V a ϕ 1 + ϕ + ϕ 2 ϕ 3 1 ϕ 3 ,
where V a = ( 4 π / 3 ) a 3 is the particle volume, and ϕ is the particle volume fraction. This expression holds to high accuracy up to the freezing transition volume fraction ϕ f 0.494 . In Figure 2b, the maximal fluid-state osmotic pressure reached at freezing, Π ( ϕ f , a ) , is plotted as a function of the particle radius a in conjunction with the thermal pressure k B T / V a , where k B is the Boltzmann constant, and T is the dispersion temperature. Notice that the 1 / a 3 dependence of both quantities implies that the osmotic pressure is strongly reduced with increasing particle size.
To illustrate that the Carnahan–Starling expression for neutral hard spheres is a decent description also for aqueous solutions of globular proteins near the isoelectric point; in Figure 3, we depict osmotic pressure data by Vilker et al. [30] as a function of ϕ , for bovine serum albumin solutions at three different pH values. As noted in [30], the isoelectric pH value is about 4.72 in a 0.15 M saline solution. We determined the volume fractions ϕ in Figure 3 from the protein concentration values and the equivalent spherical radius value a = 3.13 nm given in [30]. The pressure data at pH 4.5 are located below the Carnahan–Starling curve, suggesting that there is an attractive interaction contribution.
The gradient diffusion coefficient, D, i.e., the long-time collective diffusion coefficient, is only slightly smaller than the associated short-time collective diffusion coefficient. Hence, we can approximate D to good accuracy by its short-time form equal to the right-hand-side of [27,29,31]
D ( ϕ ) D 0 K s e d ( ϕ ) S ( ϕ ) ,
where K s e d is the short-time sedimentation coefficient, D 0 is the single-particle diffusion coefficient, and S ( ϕ ) is the osmotic compressibility factor. For the sedimentation coefficient of Brownian hard spheres, we use the analytic expression [31]
K s e d ( ϕ ) = 1 6.5464 ϕ 1 3.348 ϕ + 7.426 ϕ 2 10.034 ϕ 3 + 5.882 ϕ 4 .
As shown in Figure 4a, this expression is in excellent agreement with according dynamic simulation data where many-particles hydrodynamic interactions are accounted for. For the inverse of the osmotic compressibility factor, S ( ϕ ) , the accurate Carnahan–Starling-type expression 1 / S ( ϕ ) = ( 1 + 2 ϕ ) 2 + ϕ 3 ( ϕ 4 ) / ( 1 ϕ ) 4 is used [29]. Figure 4b displays the according (short-time) gradient diffusion coefficient, which grows monotonically with increasing particle volume fraction. The gradient diffusion coefficient, D ( ϕ ) , should not be confused with the self-diffusion coefficient, D S ( ϕ ) . In fact, the latter coefficient decreases with increasing concentration below its infinite dilution value, D 0 , both for attractive and repulsive interactions.
The dispersion shear viscosity, η , with [27]
η ( ϕ ) = η ( ϕ ) + Δ η ( ϕ ) ,
is the sum of a high-frequency (short-time) viscosity part, η , and a shear relaxation viscosity part, Δ η . The latter viscosity contribution accounts for the non-instantaneous, visco-elastic response of the dispersion microstructure to the locally generated shear flow. For Brownian hard spheres at low P e a , the two viscosity contributions are approximated, to good accuracy up to the freezing volume fraction, by [31]
η ( ϕ ) / η s = 1 + 5 2 ϕ 1 + S ^ ( ϕ ) 1 ϕ ( 1 + S ^ ( ϕ ) )
Δ η ( ϕ ) / η s = D 0 D s ( ϕ ) 12 ϕ 2 ( 1 7.085 ϕ + 20.182 ϕ 2 ) 5 ( 1 ϕ / 0.64 ) ,
with the Saito-type function S ^ ( ϕ ) = ϕ 1 + 0.95 ϕ 2.15 ϕ 2 . Here, D s ( ϕ ) / D 0 = 1 1.8315 ϕ 1 + 0.12 ϕ 0.70 ϕ 2 is a good analytic approximation for the short-time self-diffusion coefficient of Brownian hard spheres. Up to a factor of six, D s ( ϕ ) quantifies the initial slope of the mean squared displacement of hydrodynamically interacting particles. In Figure 4c, the concentration-dependent viscosity expression according to Equations (5)–(7) is depicted with simulation and experimental data. There is good agreement in the fluid-phase regime of the dispersion.
For cross-flow operating conditions compatible with the considered channel geometries, there is dispersion flow without swirling. The dispersion-averaged velocity field, V , is accordingly of the form
V ( y , z ) = v ( y , z ) y ^ + u ( y , z ) z ^ ,
where z ^ is the unit vector in axial direction. For the CM channel, y ^ is the radial unit vector pointing from the inner axis to the cylindrical membrane wall, whereas for the FMM and FMS channels, y ^ is the unit vector in the transversal direction, perpendicular to the parallel sheets (c.f. Figure 1). Furthermore, v ( y , z ) is the transversal, and u ( y , z ) is the axial velocity component.
The membrane is assumed to be fully particle retentive to the Brownian spheres. Tangential stick fluid boundary conditions are used at the inner membrane wall(s), since, in UF, the axial (i.e., z direction) solvent flow velocity inside a membrane is distinctly smaller than the axial dispersion flow velocity inside the lumen.

3. Modeling Concentration-Polarization in Ultrafiltration

Consider a dispersion of mono-disperse Brownian spheres steadily pumped through a CM, FMM, or FMS conduit as illustrated in Figure 1. We describe the dispersion flow on a coarse-grained level where the particle size and the porous structure of the membrane remain unresolved. The mass and momentum transport in UF are then described by macroscopic continuum mechanics equations governing the spatio-temporal evolution of the dispersion-averaged volume concentration field ϕ ( r , t ) of particles at position r and time t, and the dispersion-averaged velocity field V ( r , t ) .
For conditions met in UF where the single-particle Reynolds number is small compared to one, the dispersion-averaged incompressible laminar flow is described by the quasi-stationary effective Stokes equation [26,27],
P = η ( ϕ ) Δ V + d η d ϕ ϕ · V + ( V ) T
· V = 0 ,
in conjunction with the incompressibility constraint, · V = 0 , for the dispersion velocity field V . The dispersion-averaged pressure field, P ( r ) = Π ( r ) + p f ( r ) , is the sum of the equilibrium particle osmotic pressure, Π , and a fluid-phase pressure contribution, p f , adjusting itself such that incompressibility is fulfilled. There is a contribution to the Stokes equation proportional to d η / d ϕ , operative in the inhomogeneous CP layer region.
The dispersion-averaged, stationary local particle flux, J ( r ) , with
J = ϕ V D ( ϕ ) ϕ
has an advection contribution, ϕ ( r ) V ( r ) , and a diffusion contribution, D ( ϕ ( r ) ) ϕ ( r ) , respectively, with the latter quantified by the gradient diffusion coefficient D ( ϕ ) introduced in Equation (3). Substituting this flux into the macroscopic continuity equation expressing mass conservation leads to the steady-state advection–diffusion equation,
V · ϕ = · D ( ϕ ) ϕ ,
where ϕ / t = 0 . Since we assume fully particle-retentive membranes, it holds that J · y ^ = 0 at the lumen side of the membranes.
As it is discussed in [15], in the membrane interior, the pressure gradient in the axial direction is much smaller than in the transversal direction. This allows for a transversal integration of the local Darcy equation, describing the pore-size-averaged flow inside the membrane, across the membrane thickness h. The result of this integration is the Darcy–Starling (boundary condition) relation [2]
v w ( z ) = L p Δ T P ( z ) Π ( ϕ w ( z ) ) ,
for the transversal flow velocity (permeate flux), v w ( z ) , at the inner membrane wall(s). The permeate flux is here defined with a positive sign. In FMS geometry, v w is obviously zero at the impermeable lower substrate wall. Here, Δ T P ( z ) = P w ( z ) P p e r m is the local transmembrane pressure at axial distance z from the inlet, and P w ( z ) and ϕ w ( z ) are the lumen-side pressure and particle volume concentration at the membrane wall(s), respectively. The pressure at the permeate side, P p e r m , is taken as constant.
The hydraulic permeability, L p , of a clean membrane is given by [38,39]
L p = κ / η s R ln ( 1 + h / R ) ( CM ) κ / η s h ( FMM / FMS ) ,
where κ is the mean Darcy permeability of the membrane [40,41,42], averaged over its thickness h. Notice the logarithmic curvature correction for a cylindrical membrane (CM), which matters since R and h are of comparable magnitude in UF [8].
The boundary conditions at the inlet and outlet ports are, for | y | R ,
ϕ ( y , z = 0 ) = ϕ b P ( y , z = L ) = P L
and
u ( y , z = 0 ) = u 0 1 y 2 R 2 .
Here, ϕ b is the feed concentration of particles at the inlet, and P L is the pressure given at the outlet. The velocity field at the inlet port (where z = 0 ) is taken here as a fully developed parabolic flow in axial direction, so that v ( y , z = 0 ) = 0 . Furthermore, u 0 = u ( y = 0 , z = 0 ) is the axial velocity at the inlet center.
The present UF model is specified by the input (operating) parameters { ϕ b , u 0 , P L , P p e r m } . The pressure at the outlet is fixed to P L = 1 atm, i.e., it is taken as the atmospheric pressure. In the results presented in Section 5, the values of these parameters are selected such that there is no unwarranted axial flow exhaustion or permeate flow reversal. Conditions for the absence of these phenomena are discussed in [15].
For the given values of ϕ b , u 0 , and P L = 1 atm, and instead of specifying the permeate pressure, we use alternatively as a fourth input parameter the channel-length-averaged linearized transmembrane presssure, Δ T ( l ) P , referred to as linearized TMP for short. The linearized TMP is the length-average of the linear axial pressure profile minus the constant permeate pressure, i.e.,
Δ T ( l ) P ( z ) = P ( l ) ( z ) P p e r m = P ( l ) ( 0 ) P ( l ) ( 0 ) P L z L P p e r m ,
associated with a Hagen–Poiseuille-type quadratic velocity field inside the CM, FMM, and FMS channels occurring for L p = 0 and ϕ b = 0 . The brackets denote the channel-length average,
= 1 L 0 L d z .
The linearized TMP and the actual TMP, Δ T P , are thus expressed as
Δ T ( l ) P = 1 2 P ( l ) ( 0 ) + P L P p e r m
Δ T P = P ( z ) P p e r m ,
respectively. Here, P ( z ) is the, in general, non-linear dispersion pressure profile in UF. Moreover, it is P ( l ) ( 0 ) = u 0 λ 1 η s L / R H 2 + P L , with the hydraulic radius R H = { R / 2 , R , R } and the dimensionless geometry coefficient λ 1 = { 1 , 2 , 2 } for { CM , FMM , FMS } , respectively (see Table 1). The reason for using the linearized TMP as input parameter is that for values of L p commonly encountered in UF, the actual TMP is practically equal to the linearized one (see [15]). The length-averaged permeate flux is related to the TMP by v w = L p Δ T P Π .
A useful process indicator characterizing the filtration efficiency is the solvent recovery indicator, β , given by [2,43]
β = Q p e r m Q 0 = M v w A u ¯ 0 .
Here, u ¯ 0 is the cross-sectional average of the inlet velocity u ( y , z = 0 ) . The overline indicates the cross-sectional average,
¯ = 1 A A d S ,
with the constant cross-sectional area equal to A = π R 2 for CM, and A = R W both for FMM and FMS. Furthermore, Q 0 = A u ¯ 0 is the dispersion volume flow rate through the inlet cross-section, and Q p e r m = M v w is the permeate volume flow rate through the inner membrane area M. This area is equal to { 2 π R L , 2 W L , W L } for { CM , FMM , FMS } , respectively. A larger value for β reflects a larger concentration of particles at the outlet port. Notice that volume conservation implies that Q 0 = Q L + Q p e r m , where Q L is the dispersion volume flow rate through the outlet cross section.

4. Boundary Layer Analysis

In this section, we generalize the modified boundary layer analysis (mBLA) method, introduced in [15] for the CM geometry, to the flat sheets systems FMM and FMS. Since the mBLA method is explained in [15], we only summarize the essentials of the method, with our focus set on the differences between the considered flat-sheet and cylindrical membrane systems.
For standard UF operating conditions, the CP layer is a thin boundary layer of characteristic thickness δ C P R , across which ϕ is steeply decreasing, from its maximal value, ϕ w ( z ) , attained at the wall towards its minimal bulk value ϕ b . On introducing the smallness parameter ϵ δ = δ C P / R with ϵ δ R / L 1 , and for appropriately selected base units, the advection–diffusion equation is seen to be singularly perturbed. From a dominant balance analysis of this equation, the smallness parameter is identified as the inverse of the transversal Pèclet number, P e R , i.e., [15,25,44]
ϵ δ = 1 P e R = D 0 R L p Δ T ( l ) P .
The transversal Pèclet number is the ratio of the diffusion time of particles, R 2 / D 0 , and the transversal flow advection time, R / L p Δ T ( l ) P , across the transversal distance R. It is hereby assumed that the feed concentration at the inlet is small, i.e., ϕ b 1 . For a significantly developed CP layer, P e R is of the order O 10 2 . We describe in the following how the flow and concentration fields outside and inside the boundary layer are obtained and asymptotically matched.

4.1. Outer Solution

The partial differential equations determining the outer flow and concentration solutions are obtained from expanding the effective Stokes, continuity, and advection–diffusion Equations (9), (10) and (12), respectively, up to zeroth order in ϵ δ . The result is [15]
v y + u z = v y
P y = 0
P z y η ( ϕ ) u y = η ( ϕ ) y u y
v ϕ y + u ϕ z = 0 .
The curvature-related contributions on the right-hand side are non-zero for the CM geometry only. Using the associated boundary conditions specified in Section 3 up to the zeroth order in ϵ δ , and the Darcy–Starling expression, the above set of linear differential equations is solved by separation of variables and variation of constants. Depending on the channel geometry, the outer concentration and flow solutions are ϕ ( y , z ) = const . = ϕ b and
u ( y , z ) = U o u t ( y ) u 0 ( z )
v ( y , z ) = V o u t ( y ) v w ( z )
with longitudinal velocity factor
U o u t ( y ) = 1 y R 2 ( CM / FMM / FMS ) ,
and transversal velocity factor
V o u t ( y ) = 2 y R y R 3 ( CM ) 1 2 3 y R y R 3 ( FMM ) 1 4 3 y R y R 3 + 2 ( FMS ) .
To zeroth order in ϵ δ , the outer velocity field is factorized in y and z. The permeate flux, v w ( z ) , depends on the pressure field, P ( z ) , through the Darcy–Starling expression, while the axial velocity at the center-line, u 0 ( z ) = u ( 0 , z ) , is determined by the pressure according to u 0 ( z ) = λ 1 R H 2 / η s d P / d z .
The pressure, in turn, is obtained as
P ( z , [ ϕ w ] ) P p e r m Δ L ( l ) P = B + [ ϕ w ] + g ( z , [ ϕ w ] ) e K z / L + B [ ϕ w ] + g + ( z , [ ϕ w ] ) e K z / L
where Δ L ( l ) P = P ( l ) ( 0 ) P L and
g ± ( z , [ ϕ w ] ) = ± K 2 L Δ L ( l ) P 0 L e ± K z / L Π ( ϕ w ( z ) ) d z ,
B ± [ ϕ w ] = 1 2 cosh ( K ) P L P p e r m Δ L ( l ) P 1 K e K g + ( L , [ ϕ w ] ) e K g ( L , [ ϕ w ] ) e K .
We have introduced here the dimensionless effective permeability parameter, K, given by
K 2 λ 1 λ 2 η s L p L 2 R H 3 = λ 1 M L U ¯ o u t A R H · κ H R H .
It quantifies the overall solvent permeability of the membrane and the, in general, non-linear longitudinal pressure drop along the length of the membrane (c.f. [15]). The permeability parameter depends on the system geometry through the dimensionless coefficients λ 1 and λ 2 , with λ 2 = { 2 , 3 / 2 , 3 / 4 } for { CM , FMM , FMS } , on the (curvature-corrected) membrane thickness given by H = R ln 1 + h / R for CM, and H = h for FMM and FMS, respectively, and on the cross-sectional average, U ¯ o u t , of U o u t ( y ) equal to 1 / 2 for CM and 2 / 3 for FMM and FMS. Table 1 summarizes the geometry-dependent quantities determining the outer solution.
For the pressure P and transversal velocity v, no distinction is required between inner and outer solutions, since to the first order in ϵ δ , these quantities do not change steeply across the CP layer. Moreover, to the first order in ϵ δ , the pressure is independent of y.

4.2. Inner Solution

We present next the essential steps leading to the inner boundary layer solution. As noted before, by a dominant balance analysis of the advection–diffusion equation, the smallness parameter is identified as ϵ δ = 1 / P e R . The leading-order continuity, effective Stokes, and advection–diffusion equations determining the inner solutions are obtained as [15]
0 = v y
0 = y η ( ϕ ) u y
v ϕ y = y D ( ϕ ) ϕ y ,
respectively. The first (continuity) and second (Stokes) equation state that the transversal velocity, v, and the shear stress, τ = η u / y , are independent of y inside the CP layer, to leading order in ϵ δ . Thus, τ = τ w ( z ) = η ϕ w ( z ) u / y w ( z ) is the shear stress at the membrane or substrate wall. Equation (38) describes the balance of transversal advection and diffusion currents. In conjunction with the boundary conditions at the membrane and substrate walls, the inner solutions for the concentration and velocity fields are obtained as
ϕ ( y , z , [ ϕ w ] ) = ϕ w ( z ) e P e R Y ( y , z , [ ϕ w ] )
u ( y , z , [ ϕ w ] ) = τ w ( z , [ ϕ w ] ) y R 1 η ( ϕ ( y , z , [ ϕ w ] ) ) d y
v ( z , [ ϕ w ] ) = v w ( z , [ ϕ w ] ) ,
where the functional dependence on the (up to this point) unknown concentration profile, ϕ w ( z ) , is indicated. We have introduced the dimensionless function
Y ( y , z , [ ϕ w ] ) = v w ( z , [ ϕ w ] ) L p Δ T ( l ) P 1 R y R D 0 D ϕ ( y , z , [ ϕ w ] ) d y ,
appearing in the inner solutions for ϕ and u. The transversal variable y is restricted to the interval [ R , R ] for the FMS geometry, which lacks reflection symmetric with respect to the plane y = 0 . For CM and FMM, y is further restricted to [ 0 , R ] , since, for FMM, it holds that ϕ ( y ) = ϕ ( y ) as noticed in Figure 1.

4.3. Asymptotic Matching and Particle Conservation

Up to now, we determined the inner- and outer-flow solutions except for the wall concentration profile ϕ w on which they are functionally dependent. Next, the outer limit of the inner solution is asymptotically matched to the inner limit of the outer solution, using additive and multiplicative mixing rules [45] as detailed for the CM geometry in [15]. This leads to the matched asymptotic solutions
ϕ ( y , z , [ ϕ w ] ) = ϕ w ( z ) ϕ b e P e R Y ( y , z , [ ϕ w ] ) + ϕ b 1 P e R Y ( y , z , [ ϕ w ] ) e P e R Y ( y , z , [ ϕ w ] )
u ( y , z , [ ϕ w ] ) = u 0 ( z ) U ( y , z , [ ϕ w ] )
v ( y , z , [ ϕ w ] ) = v w ( z ) V o u t ( y ) .
The fields P ( z , [ ϕ w ] ) and V o u t ( y ) require no asymptotic matching and are given in Section 4.1. The longitudinal velocity factor appearing in Equation (44) for the longitudinal velocity u reads
U ( y , z , [ ϕ w ] ) = 1 + y R y R η 1 ( ϕ ( y , z , [ ϕ w ] ) ) d y 0 R η 1 ( ϕ ( y , z , [ ϕ w ] ) ) d y ,
with U ( y = 0 , z , [ ϕ w ] ) = 1 . For constant viscosity, it reduces to U ( y , z , [ ϕ w ] ) = U o u t ( y ) , where U o u t ( y ) is the parabolic profile in Equation (30).
The remaining task is to determine the particle concentration profile, ϕ w ( z ) , at the membrane wall. To this end, we employ, as a global condition, the cross-sectional particle-flux conservation law,
J z ¯ ( z , [ ϕ w ] ) = J z ¯ ( z = 0 ) J z ¯ 0 ,
where J z = J · z ^ is the longitudinal component of the particle flux in Equation (11). The overline denotes the cross-sectional average introduced in Equation (22). The above conservation law is a consequence of · J = 0 and of assuming a fully particle retentive membrane. We ignore here a small longitudinal diffusion flux contribution of O [ ϵ δ 2 ] , with the implication that J z ϕ u . Notice that J z ¯ 0 is determined, according to the inlet boundary conditions, by ϕ b and u 0 = u ( y = 0 , z = 0 ) . We deviate here from our earlier mBLA work in [15] where the inlet pressure was prescribed as an inlet boundary value instead of u 0 .
Using particle-flux conservation combined with a fixed-point iteration method adapted to the flat sheets FMM and FMS geometries, ϕ w ( z ) is numerically determined. The concentration and flow fields follow then from substituting ϕ w ( z ) into the respective matched asymptotic solutions. This constitutes our mBLA method, generalized to the FMM and FMS membrane geometries.
To analyze general features of the mBLA solution, one often introduces simplifications. Consider, for example, the longitudinal velocity factor, U ( y , z , [ ϕ w ] ) in Equation (46), entering into the matched asymptotic expression for the axial velocity u. In the bulk region away from the membrane wall(s), this factor practically equals the parabolic profile 1 y 2 / R 2 . Differences from the parabolic profile are visible in the CP boundary layer region, provided they are enhanced by dividing U ( y , z , [ ϕ w ] ) through ϵ δ . Figure 5a shows a simplified mBLA result for U ( y , z , [ ϕ w ] ) , obtained by neglecting for simplicity the z-dependence of the matched solutions in Equations (43)–(46), and using constant ϕ w ( z ) = 0.4 , ϵ δ = 1.28 × 10 2 , a = 10 nm, ϕ b = 0.001 , and R = 0.5 mm. For the suspension transport properties D ( ϕ ) and η ( ϕ ) of Brownian hard spheres, we use the accurate expressions in Equations (3) and (5). In Figure 5a, the curves of U for the symmetric geometries CM and FMM are identical and practically equal to 1 y 2 / R 2 except for the boundary layer region y ( ± ) R magnified in the inset. In the inset, U ( y , z , [ ϕ w ] ) / ϵ δ is plotted as a function of the stretched distance, ( y + R ) / δ C P = ( y / R + 1 ) / ϵ δ , from the bottom wall (in case of FMM and FMS). Regarding the FMS geometry, a slight deviation from symmetry with respect to the y = 0 plane is noticeable in the bulk region adjacent to the lower impermeable substrate wall. The inset reveals for the FMS geometry that its stretched longitudinal velocity factor near the substrate wall is larger than that for FMM, due to the absence of a CP layer at the impermeable substrate wall.
For comparison, in Figure 5b, the velocity factor V ( y ) = V o u t ( y ) = V i n ( y ) according to Equation (31) is shown for the three geometries. Note that no distinction is required between inner and outer forms of this velocity factor. Different from flat sheet geometries and owing to curvature, V ( y ) for CM has its peak value attained away from the membrane wall at radial distance y = 2 / 3 R from the center-line. While V ( y = 0 ) = 0 in CM and FMM due to symmetry, the curve of V ( y ) in FMS is strongly asymmetric with V ( y = R ) = 0 , which reflects the impermeability of the substrate wall.
This completes our presentation of the mBLA method for cylindrical and flat sheets membrane geometries. Table 2 summarizes the conditions for the validity of this method.

4.4. Remarks on the Generalized mBLA Method

Having established the generalization of the mBLA method to flat sheets geometries, a few remarks are in order here about the pros and cons of the method, its relation to other works on UF, and possible extensions to more complex dispersions.
The mBLA method for the CM geometry was shown (in [15]) to be in excellent accord with elaborate finite element calculations. While no FEM calculations are presented for the flat sheet geometries, it can be expected that the accuracy of the mBLA remains comparably good. For the CM geometry, the C++ implementation of the mBLA method is typically a thousand times faster than the according FEM calculations and it is numerically stable [15]. By using the generalized mBLA method, spatially resolved UF concentration and flux profiles are determined in addition to their spatial averages, based on theoretically precise and experimentally well-tested expressions for the dispersion properties. In case one is interested in spatial averages only, e.g., in order to compare with experimental data for the length-averaged permeate flux and TMP, a number of more simple methods are given in the literature as summarized in [2,46]. The mBLA method in its present form deals with generic CP layer effects only. Specific fouling mechanisms are not considered so far, but work on the inclusion of reversible cake layer effects is in progress. Moreover, we are currently extending the mBLA method to structured membranes as discussed in [47,48,49].
Since only the osmotic pressure and the transport properties D ( ϕ ) and η ( ϕ ) are required as input characterizing the filtered dispersion, the generalized mBLA method is readily applicable to more complex dispersions of practical relevance. These include dispersions of soft and solvent-permeable colloidal particles such as micro- and nanogels [15,31,43,50,51], and charged rigid colloidal particles, ionic microgels, and proteins [52,53,54,55]. For charged particle dispersions at lower ionic strength, the concentration dependencies of D, η , and Π are distinctly different from those of electrically neutral particles, with accordingly strong differences in the UF performance (see [50,51,54]).
In closing this section, we note that a key ingredient of the generalized mBLA method is a transversal advection–diffusion equation for the CP layer region. This equation is obtained from a dominant balance analysis explained in great detail in [15]. Such a dominant balance analysis was invoked also in earlier work by Denisov [25], but only for constant values of the gradient diffusion coefficient and dispersion viscosity and in conjunction with the highly approximate linear Van’t Hoff expression for the osmotic pressure. A dominant balance analysis of the advection–diffusion balance in the CP layer different from the one used in this work leads to another similarity solution for the CP layer profiles. The similarity solution was used in previous boundary layer theory studies of UF [56,57,58,59]. For the operating conditions employed in [15] and the present work, the similarity solution results, however, are in distinctly less good agreement with finite-element benchmark calculations (c.f. Figure 6 in [15]).

5. Results and Discussions

We quote first the employed input (operating) parameters, namely, the feed volume concentration, ϕ b , axial velocity at the inlet center, u 0 , and linearized TMP, Δ T ( l ) P . The pressure at the outlet is fixed to P L = 1 atm. These parameters are complemented by the selected membrane and dispersion parameters including R and L, membrane hydraulic permeability, L p , solvent viscosity, η s , of water at room temperature, and particle radius a. While these parameters are partially varied in the presented results, as a reference set of values, we use ϕ b = 0.001 , R = 0.5 mm, L = 0.5 m, and L p = 6.7 × 10 10 m/(Pa s). The selected value for L p is in the range of typical UF permeabilities [2], and a = 3.13 nm. Reference values for u 0 depend on the membrane geometry and are selected as { 6.80 , 5.09 , 5.09 } × 10 2 m/s for { CM , FMM , FMS } , respectively. In all three geometries, the same value, u ¯ 0 = 3.40 × 10 2 m/s, is obtained for the cross-section averaged inlet velocity, where u CM 0 = 4 / 3 × u FMM / FMS 0 and u 0 = u ( 0 , 0 ) . The linearized TMP is varied from 1 kPa up to 20 kPa. The filtrated dispersion consists of mono-disperse Brownian hard spheres, suspended in water at room temperature. The particle transport properties, D ( ϕ ) and η ( ϕ ) , and the osmotic pressure, Π ( ϕ ) , are accounted for in the mBLA method using accurate analytic expressions given in Section 2.
Using the reference values for { R , L , L p , η s } , the geometry-dependent reference values for the effective permeability parameter K defined in Equation (35) are obtained as { 0.146 , 0.0634 , 0.0448 } for { CM , FMM , FMS } , respectively. Since K 2 1 , the actual TMP, Δ T P , is practically equal to the linearized TMP = Δ T ( l ) P . For this reason, we do not distinguish any more in the following between actual and linearized TMP.
In Section 5.1, we discuss mBLA results for the wall concentration profile, ϕ w ( z ) , and for the longitudinal particle flux J z ( y , z ) . For the latter quantity, both its bulk and excess parts are discussed. The TMP is not varied here, but two different particle radii a = 3.13 and 10 nm are considered. In Section 5.2, the TMP, feed concentration, and cross-section averaged inlet velocity dependence of the length-averaged permeate flux, v w , and of the solvent recovery indicator β defined in Equation (21) are studied. Triggered by the results in Section 5.2, in Section 5.3, three independent parameter combinations out of the vast number of operating, membrane, and dispersion parameters are identified, which solely determine the solvent recovery indicator β .

5.1. CP-Layer and Longitudinal Particle Transport for Reference Conditions

For the UF geometries in Figure 6, we compare the concentration profile, ϕ w ( z ) , at the membrane wall for particle radii (a) a = 3.13 nm and (b) a = 10 nm, respectively, for fixed transversal Pèclet number P e R = 78 , and for the reference values of L p , R, and L. Consider first Figure 6a for the concentration profile, where the reference value for the average inlet velocity, u ¯ 0 = 3.40 × 10 2 m/s, is used in all three geometries, in addition to a common TMP value of 16 kPa. The shape of the concentration profile is quite similar in the three cases, with the maximal concentration values at the outlet being distinctly smaller than the freezing transition volume concentration ϕ f = 0.494 of colloidal hard spheres. The concentration profiles for FMM and FMS are nearly the same along the full channel length. Slight differences between these profiles, and the crossover of the CM with the other two profiles at axial distance z / L > 0.6 , can be qualitatively understood from inspecting (in the inset) the axial velocity, u ( y = 0 , z ) , at the axial center-line y = 0 . Note that for the CM geometry u 0 is larger than for FMM and FMS, by a factor of 4 / 3 . At axial distances similar to those for the ϕ w ( z ) profile, crossovers are observed between the three curves for u ( 0 , z ) . Since d u ( 0 , z ) / d z = ( λ 2 / R H ) v w ( z ) , a stronger decay of u ( 0 , z ) reflects a larger transmembrane flux, v w ( z ) , and thus an enhanced transport of particles towards the membrane wall.
The usage of the same value for u ¯ 0 in the three geometries implies that the shear rate at the inlet of the membrane wall, γ ˙ * = ( d u ( y , z = 0 ) / d y ) ( y 0 ) = 2 u 0 / R (c.f. Equation (1)), is larger for CM than for FMM and FMS. This explains the initially smaller increase in ϕ w ( z ) for CM. The dashed (red) curve in Figure 6a holds when for FMM the same value for the wall shear rate γ ˙ * is used as for CM (FMM2 curve). While the initial slope of ϕ w ( z ) at the inlet is the same for CM and FMM2, the ϕ w ( z ) values for FMM2 sufficiently distant from the inlet are distinctly smaller than in the CM case. This is reflected in the axial velocity profile u ( 0 , z ) depicted in red in the inset.
A similar analysis as for Figure 6a applies to Figure 6b, for a dispersion of larger particles of radius a = 10 nm, smaller TMP = 5 kPa, and cross-section averaged inlet velocity u ¯ 0 = 1.11 × 10 2 m/s equal to that of CM, FMM, and FMS. The value for γ ˙ * in FMM2 is equal to the one in CM. While the same general trends are encountered as in Figure 6a, there are now distinctly larger wall concentrations observed, with ϕ w ( L ) > 0.32 at the channel outlet.
The ratio, Π / TMP , of length-averaged transmembrane osmotic pressure and TMP equals 0.34 and 0.18 in Figure 6a,b, respectively. It reflects the strong influence of the CP layer on the permeate flux v w , exerted through the (length-averaged) Darcy–Starling law in Equation (13). Notice further that for extreme parameter values not considered here where the longitudinal pressure drop across the channel length is comparable to the TMP, the according profiles for ϕ w ( z ) can be non-monotonic in z, with maximum concentration attained for z < L . As discussed in [15], in the context of the CM geometry, a criterion for a strictly monotonic increase in ϕ w ( z ) is that TMP > 3 P ( 0 ) P L .
We proceed by discussing the axial particle flux, J z ( y , z ) , which according to
J z ( y , z ) = ( ϕ ( y , z ) ϕ b ) u ( y , z ) + ϕ b u ( y , z ) J z e x ( y , z ) + J z b ( y , z ) ,
consists of an excess flux contribution, J z e x ( y , z ) , that vanishes in the bulk region and a bulk flux contribution, J z b ( y , z ) . Right at the membrane wall, J z e x ( R , z ) = 0 due to the employed stick boundary condition. Without a significant CP layer, J z e x is practically zero. The results by the mBLA method for the excess and bulk parts of J z ( y , z ) are depicted in Figure 7 for CM and FMM, respectively, using a = 3.13 nm and input parameters as in Figure 6a. The fluxes are plotted semi-logarithmically as functions of the (relative) distance, 1 y / R , from the (upper, in case of FMM) membrane wall. The FMS fluxes are not shown here since near the upper membrane, they are practically equal to those of the FMM. The excess flux is zero at the impermeable substrate wall of the FMS channel. The displayed flux profiles are given in units of the cross-sectional flux average, J z ¯ 0 = ϕ b u ¯ 0 , at the inlet, which according to Equation (47), is z-independent. The excess axial flux (solid curves in Figure 7) increases from zero at the membrane ( y = R ) towards its maximal value at y δ C P . The characteristic CP layer thickness, δ C P = R / P e R , is marked by the vertical dotted lines. Due to a trade-off between ϕ and u, where the former decreases and the latter increases with increasing distance from the membrane, the maximum of J z e x is located at a distance from the wall larger than δ C P , in the transition region between CP layer and bulk. The flux maximum grows with increasing distance from the inlet, and its location shifts away from the membrane wall. There is an according decrease, with increasing z, of the maximum values of J z b at the channel center-line y = 0 .
According to Figure 7, the excess fluxes in CM are similar to those in FMM, which can be attributed to the fact that curvature effects in the thin CP layer of the CM geometry are negligible. In contrast, there are pronounced differences observed for the maximal bulk fluxes of J z b ( y = 0 , z ) = ϕ b u ( 0 , z ) , reflected in the likewise different behavior of u ( 0 , z ) for CM and FMM, seen in the inset of Figure 6a.
To elucidate differences in the axial particles transport between CM and the flat sheets FMM and FMS geometries, in Figure 8, the cross-section averaged axial excess flux, J z ¯ e x ( z ) , is shown for a = 3.13 nm and a = 10 nm, respectively. Like ϕ w ( z ) , the average flux increases monotonically from 0 at the inlet towards its maximal value at the outlet where the CP layer becomes most pronounced. In both panels of Figure 8, where the same value for u ¯ 0 is used, the excess flux for CM is twice as large as the one for FMM. To understand this, notice that the cross-section averaged axial bulk flux follows, from the sum rule expressing particle flux conservation and our matched asymptotic solution, as
J z ¯ b ( z ) / J z ¯ 0 = 1 J z ¯ e x ( z ) / J z ¯ 0 = M A u ¯ 0 1 L 0 z d z v w z , [ ϕ w ] + O ϵ δ .
Recall next from Table 1 that the ratio, M / A , of membrane area M and cross-section area A is twice as large for CM as for FMM. One further observes in Figure 8 that J z ¯ e x ( z ) F M M 2 × J z ¯ e x ( z ) F M S , which arises since the value for M / A in FMS is half that of FMM, while v w ( z ) is practically the same in both cases. In comparing Figure 8a,b, one notices, for operating conditions as those in Figure 6a,b, that the normalized excess flux curves in (b) are larger than the according ones in (a). In general, there is an intermingled influence of D, η , and Π on the UF behavior when the particle size is varied.

5.2. TMP, Feed Concentration, and Velocity Effects on Global Indicators

Having discussed the local profiles ϕ w ( z ) and J z ( y , z ) , we assess next the influence of varying TMP and feed concentration on properties that globally indicate the UF performance. Two global indicators of particular importance are the length-averaged permeate flux, v w = L p TMP Π , and the solvent recovery indicator β 1 introduced in Equation (21). The latter indicator gives the fraction of initial dispersion volume, recovered as pure solvent in the permeate compartment. It is related to another indicator, α = ϕ ¯ L / ϕ b 1 , defined as the ratio of cross-section averaged concentration at the outlet and the (cross-section averaged) feed concentration ϕ b . Owing to particle flux and dispersion volume conservation valid for an ideally particle retentive membrane, the concentration factor α is related to β simply by [43]
α = 1 1 β .
One generally observes that v w increases monotonically with increasing TMP. Considering Equations (21) and (50), this implies that both β and α are monotonically increasing with increasing TMP.
Our central aim is to identify dimensionless combinations of input parameters characterizing uniquely the global process indicators, independent of the considered geometries. As a prerequisite, in Figure 9, we analyze the process indicators v w , α , and β as functions of TMP, for input parameters selected otherwise as in Figure 6a,b. The symbols are predictions by the mBLA method. Regarding the average permeate flux depicted in panels (a-1) and (b-1), the mBLA results for v w are nearly coincidental, even for the largest TMP considered. The straight solid lines depict the pure-solvent form v w 0 = L p TMP , for the same hydraulic permeability L p used in both figure panels. The influence of the CP layer (osmotic pressure) buildup on the permeate flux is noticeable at larger TMP where v w increased sub-linearly. For particle radius 10 nm, a smaller TMP interval is depicted in the figure since cake formation by freezing is observed already for TMP 7 kPa.
The straight lines in panels (a-2) and (b-2) of Figure 9 are the pure-solvent prediction, β 0 , equal to β at zero feed concentration ϕ b = 0 . Explicitly,
β 0 = M A L p u ¯ 0 TMP ,
where M is the lumen-side membrane surface area, and A is the area of the inner cross-section. As noticed from Equation (51), β 0 has a geometry-dependent factor, M / A , of values listed in Table 1, and it includes the cross-section averaged inlet velocity u ¯ 0 , which has different values in panels (a) and (b), respectively. At small TMP, the mBLA curves for β overlap with the pure-solvent lines of β 0 , since the osmotic pressure contribution is here negligible. For comparison, panels (a-3) and (b-3) show the associated concentration factor α . Solid lines represent α 0 = 1 / ( 1 β 0 ) , which quantifies the initial increase in α with increasing TMP. Notice that, according to Equation (51), both β 0 and α 0 are solely dependent on given input parameters. According to Table 1, it holds that β CM 0 : β FMM 0 : β FMS 0 = 4:2:1 for equal TMP.
The effect of a varying radius a on the length-averaged reduced permeate flux and on α is analyzed in Figure 10 (left) and (right), respectively. The permeate flux is rendered non-dimensional by multiplication with R / D 0 , and it is plotted versus ϕ b multiplied by P e R 78 , i.e., the selected value for the transversal Pèclet number. The multiplication of ϕ b by P e R is used here to highlight that ϕ b P e R < 1 , an inequality required for the validity of the mBLA method (c.f. Table 2). Except for ϕ b , all input parameters are as in Figure 6a,b, respectively. Regarding Figure 10 (left), one notices that v w ( R / D 0 ) P e R for ϕ b 0 . According to the inset, the initial plateau value of the reduced flux, equal to P e R , is mirrored by the plateau value Π / TMP 0 . With further increasing ϕ b , the accordingly increasing Π causes a decline in the reduced permeate flux linear in the logarithm of ϕ b , which was also observed in typical UF experiments [9]. For fixed a, the reduced flux is nearly the same for the considered geometries, except when ϕ b P e R is close to one. The mBLA results for α as a function of ϕ b are displayed in Figure 10 (right). The concentration factor decreases with increasing ϕ b , and it is distinctly larger for the CM geometry. One notices further in Figure 9(a-1,b-1) that different channel geometries have little influence on the permeate flux. Figure 10 (right) is useful for re-circulation setups where the concentrated dispersion, collected at the outlet port, is fed back at the inlet port.
The dependence of β on the cross-section averaged axial inlet velocity, u ¯ 0 , is shown in Figure 11, for the three geometries and for two particle sizes. The velocity is made non-dimensional by multiplication with the axial velocity unit L D 0 / R 2 . With increasingly large inlet velocity, the mBLA data for β (symbols) converge towards the pure solvent prediction β 0 (lines). This convergence arises since the influence of the CP layer ceases with increasing axial flow velocity.

5.3. Universal Behavior of Global UF Indicators

In Section 5.2, we analyzed the global UF indicators α , β , and v w ( R / D 0 ) , as functions of TMP, ϕ b , a, and u ¯ 0 , respectively, with the other input parameters kept constant. It was shown that the indicators are significantly affected by changes in these parameters and in the membrane. The indicators depend in addition on R, L, and L p , which characterize the membrane geometry and hydraulic permeability, respectively.
Since the UF process depends on a vast number of input and membrane geometry parameters, it is important to identify a minimal set of dimensionless combinations of input parameters fully characterizing the UF indicators. By extensive parameter variation studies based on the mBLA method combining numerical efficiency with high accuracy, and by using a simplified mBLA solution described in Appendix A, we succeeded in identifying three independent dimensionless combinations of input parameters, termed variables in the following, on which β is solely dependent. The first of these independent variables is the transversal Pèclet number, P e R , introduced in Equation (23) as the ratio of diffusion and flow advection times for the particle transported across the transversal distance R. The remaining two variables are identified as
β 0 P e R = M D 0 A R u ¯ 0
γ 0 P e R = 9 c 2 ϕ b L p η s R a 2 β 0 P e R .
Here, β 0 is the pure solvent recovery indicator introduced in Equation (51) where its linear dependence on TMP was emphasized. While β 0 is independent of dispersion properties, these are accounted for by γ 0 and P e R . The third variable, γ 0 , includes a geometry-dependent factor c = { 1 / 8 , 1 / 3 , 2 / 3 } for { CM , FMM , FMS } , respectively. Moreover, it depends on the feed concentration ϕ b , particle radius a, hydraulic membrane permeability L p , solvent viscosity η s , and the inner channel half-height R. The appearance and the form of γ 0 are motivated by a simplified mBLA calculation of v w R / D 0 presented in Appendix A.
On expressing β in terms of the three variables, we find that
β F P e R , β 0 P e R , γ 0 P e R
is valid to good accuracy in the assessed parameter space region. Here, F is a function depending on the three variables only.
The universal behavior of β in terms of the three variables is exemplified in Figure 12, where it is plotted as a function of P e R , for values of the duplet β 0 / P e R , γ 0 / P e R given by the two sets G 1 = 8.08 × 10 3 , 1.55 × 10 4 (open symbols) and G 2 = ( 8.08 × 10 3 , 3.10 × 10 4 ) (filled symbols), respectively. Within each set, the mBLA data for β coincide practically for all P e R and for the three geometries CM, FMM, and FMS. The two curves F ( P e R ; G 1 ) and F ( P e R ; G 2 ) traced out by the mBLA data for β merge at small P e R with the straight dotted line for β 0 determined by Equation (52). The two sets G1 and G2 for β 0 / P e R , γ 0 / P e R share the same value 8.08 × 10 3 for β 0 / P e R . The input parameters for the curves G1 and G2 in Figure 12 are listed in the rows G1 and G2 of Table 3. Notice from Figure 12 that, for given (non-small) P e R and β 0 , a larger value of γ 0 causes a smaller values of β due to the lowered permeate flux. According to Equation (53), γ 0 increases for increasing ϕ b and L p , and for decreasing particle radius a.
In rows G1 and G2 of Table 3, the membrane parameters R, L, and L p are kept constant. To illustrate that changes in these parameters are likewise compatible with Equation (54), in Figure 13, β is shown as a function of P e R , for three sets of the duplet β 0 / P e R , γ 0 / P e R given by S - CM = 8.08 × 10 3 , 1.55 × 10 4 (circles), S - FMM = 4.04 × 10 3 , 2.07 × 10 4 (squares), and S - FMS = 2.02 × 10 3 , 2.07 × 10 4 (triangles). These three sets are not taken out of the blue but are precisely the values for β 0 / P e R and γ 0 / P e R of the reference system introduced in Section 5 for CM, FMM, and FMS, respectively. For each geometry with its according set for β 0 / P e R , γ 0 / P e R , four different lists of input parameters are given in Table 3, e.g., rows CM 1-4 for the CM geometry. Rows CM 1, FMM 1, and FMS 1 include the reference system input data for CM, FMM, and FMS, respectively. Different values for the membrane properties are used for the three geometries. In lists CM 2-4, e.g., the thickness h of the membrane is specified instead of L p , by assuming the same mean Darcy permeability, κ , as for the reference system with CM geometry. The hydraulic permeability follows then from Equation (14). According to Figure 13, β traces out three master curves for S-CM (top), S-FMM (middle), and S-FMS (bottom), respectively, where the according duplet values for ( β 0 / P e R , γ 0 / P e R ) are noted in the figure. The solid curves for β are obtained using reference system data listed in rows CM 1, FMM 1, and FMS 1 of Table 3. The dashed straight lines are the pure solvent predictions β 0 for β . It holds that β S - CM 0 > β S - FMM 0 > β S - FMS 0 (see also Figure 9(a-2)). The mBLA data for β in Figure 13 nicely reconfirm the validity of Equation (54).
The universal dependence of β on the three variables P e R , β 0 / P e R , and γ 0 / P e R in Equation (54) is inherited by the concentration factor, α , and by the reduced permeate flux, v w ( R / D 0 ) , since these are related to α by
β = 1 1 α = v w R D 0 β 0 P e R .
The independent variables on which β , α , and the reduced permeate flux are solely dependent are identified using a simplified mBLA calculation where the gradient diffusion coefficient and viscosity are approximated by their infinite dilution values, and the osmotic pressure is approximated by the linear van’t Hoff expression. Details of this calculation are given in Appendix A. The result of the Appendix is the simplified mBLA expression
v w R D 0 β 0 P e R + γ 0 P e R g 1 .
Here, g is a complicated function given in Equation (A6) of the Appendix A, which depends on P e R and other input parameters, and implicitly also on the reduced permeate flux. Thus, a fixed-point iteration solver is required also for the simplified mBLA expression in Equation (56). Note that Equation (56) holds for larger values of P e R only, and different from the full mBLA solution, it violates to some extent the sole dependence on β on the three variables in Equation (54). The virtue of Equation (56) is that it suggests that P e R , β 0 / P e R , and γ 0 / P e R are the reduced set of independent variables and that it provides, in addition, the expression for γ 0 in Equation (53).

6. Conclusions

We generalized the mBLA method for predicting local UF concentration and flow profiles, and global process indicators, from the cylindrical membrane geometry to flat sheets geometries with two and one membrane sheet, respectively. The semi-analytic form of the method makes it well suited for identifying general trends and unifying features, which are tasks that require extensive variations of different input parameters. Importantly, the semi-analytic mBLA method gives physical insight into the functional form of concentration and flow profiles and of global process indicators.
For simplicity, we considered dispersions of mono-disperse Brownian hard spheres, using accurate analytic expressions for the transport properties and osmotic pressure, as input to the generalized mBLA method. We elaborated on the effects of the three geometries on local concentration and particle flux profiles, by revealing significant differences in the average particle fluxes. For varying operating conditions, we analyzed the influence of the different geometries on global UF efficiency indicators including β , concentration factor α , and length-averaged permeate flux v w . We showed that the influence of different geometries is particularly pronounced for these indicators.
Our key result is that the intricate and intermingled dependencies of the UF indicators on a vast number of input parameters and the membrane geometry can be grouped together into only three independent variables: P e R , β 0 / P e R , and γ 0 / P e R . The semi-analytic form and numerical efficiency of the mBLA method allowed us to identify these dimensionless variables and to show that β is uniquely determined by them. The simple three-variables dependency carries over to related global indicators, and it can be of help in the design, optimization, and scale-up of UF setups.
We expect that a description of β along the lines of Equation (54) holds also for more complex dispersions such as solvent-permeable hard spheres and non-ionic nanogels, based on our experience gained for these systems [15,31,43]. It is likely that additional independent variables (parameter groups) come into play for these systems on which β is depending on. The generalized mBLA method can be used to identify these additional variables. In dealing with UF and particle size polydispersity, partially retentive membranes and fouling mechanisms are also of importance. These features were not accounted for in the present work focusing on generic CP layer and geometry effects. Work on extensions of the mBLA method to account for cake-layer effects and polydispersity is in progress.

Author Contributions

Conceptualization, G.W.P.; methodology, G.W.P. and G.N.; writing—original draft preparation, G.W.P.; writing—review and editing, G.W.P. and G.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Deutsche Forschungsgemeinschaft (SFB-985, Project B6) grant number 191948804.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data supporting the findings of this study are available from the corresponding author upon reasonable request. Our Python code for the mBLA method is freely avaliable at [60].

Acknowledgments

The authors gratefully acknowledge M. Brito (University of Stuttgart) and J.K.G. Dhont (Forschungszentrum Jülich) for helpful discussions.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

CPconcentration-polarization
CMcylindrical membrane
FMMflat sheet membranes (top and bottom)
FMSflat sheet membrane (top)/substrate (bottom)
mBLAmodified boundary layer analysis (method)
TMPtransmembrane pressure
UFultrafiltration

List of Symbols (SI Units)

P e a single-particle shear-Pèclet number
P e R transversal Pèclet number
α concentration factor
β solvent-recovery indicator
α 0 , β 0 the infinite dilution value of α , β
γ 0 third dimensionless variable characterizing β
ϵ δ = 1 / P e R perturbation parameter
δ C P characteristic thickness of CP layer (m)
R, L, Whalf-height, axial length, width of channel (m)
R H hydraulic radius (m)
h, Hmembrane thickness, curvature-corrected thickness (m)
A, Mchannel cross-section, membrane surface area (m2)
y, ztransversal, longitudinal coordinate (m)
V dispersion-averaged velocity (m/s)
v, utransversal, longitudinal velocity (m/s)
u 0 axial velocity at center of inlet (m/s)
U o u t longitudinal velocity factor of outer solution
Uasymptotically matched longitudinal velocity factor
v w permeate flux (m/s)
V o u t transversal velocity factor
γ ˙ * shear rate at inlet of membrane wall (1/s)
τ w shear stress at membrane wall (Pa)
Pdispersion-averaged pressure (Pa)
P p e r m , P L pressure at permeate side, at outlet port (Pa)
Δ T P length-averaged transmembrane pressure (Pa)
Δ T ( l ) P length-averaged, linearized transmembrane pressure (Pa)
Π particles osmotic pressure (Pa)
ρ dispersion mass density (kg/m3)
nparticle number density (1/m3)
aradius of hard spheres (m)
V a particle volume (m3)
ϕ particle volume fraction
ϕ b , ϕ w particle volume fraction at inlet (feed), at membrane wall
η , η s suspension-, solvent viscosity (Pa s)
Dgradient diffusion coefficient (m2)/s)
D 0 Stokes–Einstein diffusion coefficient (m2)/s)
κ mean Darcy permeability (m2)
L p hydraulic permeability of clean membrane (m/(Pa s))
Kdimensionless effective permeability parameter
Q 0 , Q p e r m , Q L volume flow rate through inlet, membrane, outlet (m3)/s)
J z , J z e x , J z b longitudinal particle-flux, excess part, bulk part (m/s)
length-average of (c.f. Equation (18))
¯ cross-section average of (c.f. Equation (22))

Appendix A. Simplified mBLA Method

We present here the derivation of the simplified mBLA expression for the reduced length-averaged permeate flux, v w R / D 0 , given in Equation (56) of Section 5.3.
On approximating the gradient diffusion coefficient, D ( ϕ ) , and the dispersion viscosity, η ( ϕ ) , by their infinite dilution values D 0 and η s , respectively, the mBLA expression for the concentration profile at the inner membrane wall, ϕ w ( z ) , is reduced to
ϕ w ( z ) ϕ b c ϕ b P e R 2 β 0 v w ( z , [ ϕ w ] ) / v * 2 u 0 ( z , [ ϕ w ] ) / u 0 1 L 0 z v w ( z , [ ϕ w ] ) / v * d z .
Here, v * = L P Δ T P , and c is a geometry-dependent factor equal to { 1 / 8 , 1 / 3 , 2 / 3 } for the CM, FMM, and FMS geometries, respectively. The functional dependence of v w ( z ) and u 0 ( z ) = u 0 ( y = 0 , z ) on the wall concentration profile is indicated by the bracket [ ϕ w ] .
Setting z = L and assuming ϕ w ϕ b , it follows in conjunction with the matched asymptotic solution relation u 0 ( L ) / u 0 = 1 β 0 v w / v * inferred from Equations (44) and (45) that
v w R D 0 P e R β 0 1 + c ϕ b P e R 2 v w ( L ) / v * 2 ϕ w ( L ) 1 P e R β 0 + γ .
The function γ introduced on the right-hand-side of the second equality is defined in terms of the left-hand-side. In an ensuing step, we invoke the linear van’t Hoff osmotic pressure expression, Π = ( k B T / V a ) ϕ , for the wall concentration ϕ w ( L ) at the outlet. This leads to
ϕ w ( L ) P e R 2 9 a 2 L P R η Π ( ϕ w ( L ) ) Δ T P .
Substitution of the above expression for ϕ w ( L ) / P e R into Equation (A2), and solving for γ , results into
γ β 0 9 2 c ϕ b L p η s R a 2 × P e R v w ( L ) v * 2 Δ T P Π ( ϕ w ( L ) ) = γ 0 g ,
where γ 0 and the function g are identified as
γ 0 = β 0 9 2 c ϕ b L p η s R a 2
g = P e R v w ( L ) v * 2 × Δ T P Π ( ϕ w ( L ) ) .
In combination with Equation (A2), this leads to Equation (56), with the function g given in Equation (A6).

References

  1. Mason, E.A.; Lonsdale, H.K. Statistical-mechanical theory of membrane transport. J. Membr. Sci. 1990, 51, 1–81. [Google Scholar] [CrossRef]
  2. Mulder, M. Basic Principles of Membrane Technology; Springer: Dordrecht, The Netherlands, 2013. [Google Scholar]
  3. Bacchin, P.; Si-Hassen, D.; Starov, V.; Clifton, M.J.; Aimar, P. A unifying model for concentration polarization, gel-layer formation and particle deposition in cross-flow membrane filtration of colloidal suspensions. Chem. Eng. Sci. 2002, 57, 77–91. [Google Scholar] [CrossRef] [Green Version]
  4. Blatt, W.F.; Dravid, A.; Michaels, A.S.; Nelsen, L. Solute polarization and cake formation in membrane ultrafiltration: Causes, consequences, and control techniques. In Membrane Science and Technology; Flinn, J.E., Ed.; Springer: Boston, MA, USA, 1970; pp. 47–97. [Google Scholar]
  5. Belfort, G.; Davis, R.H.; Zydney, A.L. The behavior of suspensions and macromolecular solutions in crossflow microfiltration. J. Membr. Sci. 1994, 96, 1–58. [Google Scholar] [CrossRef]
  6. Bacchin, P.; Aimar, P.; Field, R.W. Critical and sustainable fluxes: Theory, experiments and applications. J. Membr. Sci. 2006, 281, 42–69. [Google Scholar] [CrossRef] [Green Version]
  7. Bacchin, P.; Marty, A.; Duru, P.; Meireles, M.; Aimar, P. Colloidal surface interactions and membrane fouling: Investigations at pore scale. Adv. Colloid Interface Sci. 2011, 164, 2–11. [Google Scholar] [CrossRef]
  8. Doran, P.M. Unit operations. In Bioprocess Engineering Principles; Elsevier: Amsterdam, The Netherlands, 2003; pp. 445–595. [Google Scholar]
  9. Baker, R.W. Membrane Technology and Applications; John Wiley & Sons: Hoboken, NJ, USA, 2004. [Google Scholar]
  10. Bopape, M.; Van Geel, T.; Dutta, A.; Van der Bruggen, B.; Onyango, M. Numerical modelling assisted design of a compact ultrafiltration (UF) flat sheet membrane module. Membranes 2021, 11, 54. [Google Scholar] [CrossRef]
  11. Rey, C.; Hengl, N.; Baup, S.; Karrouch, M.; Dufresne, A.; Djeridi, H.; Dattani, R.; Pignon, F. Velocity, stress and concentration fields revealed by micro-PIV and SAXS within concentration polarization layers during cross-flow ultrafiltration of colloidal Laponite clay suspensions. J. Membr. Sci. 2019, 578, 69–84. [Google Scholar] [CrossRef]
  12. Meng, B.; Li, X. In situ visualization of concentration polarization during membrane ultrafiltration using microscopic laser-induced fluorescence. Environ. Sci. Technol. 2019, 53, 2660–2669. [Google Scholar] [CrossRef]
  13. Lüken, A.; Linkhorst, J.; Fröhlingsdorf, R.; Lippert, L.; Rommel, D.; Laporte, L.D.; Wessling, M. Unravelling colloid filter cake motions in membrane cleaning procedures. Sci. Rep. 2020, 10, 20043. [Google Scholar] [CrossRef]
  14. Pipich, V.; Starc, T.; Buitenhuis, J.; Kasher, R.; Petry, W.; Oren, Y.; Schwahn, D. Silica fouling in reverse osmosis systems-operando small-angle neutron scattering studies. Membranes 2021, 11, 413. [Google Scholar] [CrossRef] [PubMed]
  15. Park, G.W.; Nägele, G. Modeling cross-flow ultrafiltration of permeable particle dispersions. J. Chem. Phys. 2020, 153, 204110. [Google Scholar] [CrossRef]
  16. Marshall, E.; Trowbridge, E. Flow of a Newtonian fluid through a permeable tube: The application to the proximal renal tubule. Bull. Math. Biol. 1974, 36, 457–476. [Google Scholar] [CrossRef]
  17. Pozrikidis, C. Stokes flow through a permeable tube. Arch. Appl. Mech. 2010, 80, 323–333. [Google Scholar] [CrossRef]
  18. Tilton, N.; Martinand, D.; Serre, E.; Lueptow, R.M. Incorporating Darcy’s law for pure solvent flow through porous tubes: Asymptotic solution and numerical simulations. AIChE J. 2012, 58, 2030–2044. [Google Scholar] [CrossRef]
  19. Porter, M.C. Concentration polarization with membrane ultrafiltration. Ind. Eng. Chem. Prod. Res. Dev. 1972, 11, 234–248. [Google Scholar] [CrossRef]
  20. Trettin, D.R.; Doshi, M.R. Limiting flux in ultrafiltration of macromolecular solutions. Chem. Eng. Commun. 1980, 4, 507–522. [Google Scholar] [CrossRef]
  21. Berman, A.S. Laminar flow in channels with porous walls. J. Appl. Phys. 1953, 24, 1232–1235. [Google Scholar] [CrossRef]
  22. Brian, P.L.T. Concentration polarization in reverse osmosis desalination with variable flux and incomplete salt rejection. Ind. Eng. Chem. Fundam. 1965, 4, 439–445. [Google Scholar] [CrossRef]
  23. Bouchard, C.R.; Carreau, P.J.; Matsuura, T.; Sourirajan, S. Modeling of ultrafiltration: Predictions of concentration polarization effects. J. Membr. Sci. 1994, 97, 215–229. [Google Scholar] [CrossRef]
  24. Kleinstreuer, C.; Paller, M.S. Laminar dilute suspension flows in plate-and-frame ultrafiltration units. AIChE J. 1983, 29, 529–533. [Google Scholar] [CrossRef]
  25. Denisov, G. Theory of concentration polarization in cross-flow ultrafiltration: Gel-layer model and osmotic-pressure model. J. Membr. Sci. 1994, 91, 173–187. [Google Scholar] [CrossRef]
  26. Dhont, J.K.G. An Introduction to Dynamics of Colloids; Elsevier: Amsterdam, The Netherlands, 1996. [Google Scholar]
  27. Nägele, G.; Dhont, J.K.G.; Voigtmann, T. Theory of colloidal suspension structure, dynamics, and rheology. In Theory and Applications of Colloidal Suspension Rheology; Cambridge University Press: Cambridge, UK, 2021; pp. 44–119. [Google Scholar]
  28. Hansen, J.P.; McDonald, I.R. Theory of Simple Liquids; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
  29. Banchio, A.J.; Nägele, G. Short-time transport properties in dense suspensions: From neutral to charge-stabilized colloidal spheres. J. Chem. Phys. 2008, 128, 104903–104920. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Vilker, V.L.; Colton, C.K.; Smith, K.A. The osmotic pressure of concentrated protein solutions: Effect of concentration and pH in saline solutions of bovine serum albumin. J. Colloid Interface Sci. 1981, 79, 548–566. [Google Scholar] [CrossRef]
  31. Riest, J.; Eckert, T.; Richtering, W.; Nägele, G. Dynamics of suspensions of hydrodynamically structured particles: Analytic theory and applications to experiments. Soft Matter 2015, 11, 2821–2843. [Google Scholar] [CrossRef] [Green Version]
  32. Ladd, A.J.C. Hydrodynamic transport coefficients of random dispersions of hard spheres. J. Chem. Phys. 1990, 93, 3484–3494. [Google Scholar] [CrossRef]
  33. Segre, P.; Meeker, S.; Pusey, P.; Poon, W. Viscosity and structural relaxation in suspensions of hard-sphere colloids. Phys. Rev. Lett. 1995, 75, 958. [Google Scholar] [CrossRef]
  34. Abade, G.; Cichocki, B.; Ekiel-Jezewska, M.; Nägele, G.; Wajnryb, E. Diffusion, sedimentation, and rheology of concentrated suspensions of core-shell particles. J. Chem. Phys. 2012, 136, 104902. [Google Scholar] [CrossRef] [Green Version]
  35. Foss, D.R.; Brady, J.F. Structure, diffusion and rheology of Brownian suspensions by Stokesian dynamics simulation. J. Fluid Mech. 2000, 407, 167–200. [Google Scholar] [CrossRef] [Green Version]
  36. Phung, T.N. Behavior of Concentrated Colloidal Suspensions by Stokesian Dynamics Simulation. Ph.D. Thesis, California Institute of Technology, Pasadena, CA, USA, 1993. [Google Scholar]
  37. Weiss, A.; Dingenouts, N.; Ballauff, M.; Senff, H.; Richtering, W. Comparison of the effective radius of sterically stabilized latex particles determined by small-angle X-ray scattering and by zero shear viscosity. Langmuir 1998, 14, 5083–5087. [Google Scholar] [CrossRef]
  38. Stevenson, J.F.; Parry, J.S.; Gupta, K.M. Hydraulic permeability of hollow-fiber membranes. J. Biomed. Mater. Res. 1978, 12, 401–419. [Google Scholar] [CrossRef]
  39. Bird, R.B.; Stewart, W.E.; Lightfoot, E.N.; Anderson, W. Transport Phenomena; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2002. [Google Scholar]
  40. Renard, P.; De Marsily, G. Calculating equivalent permeability: A review. Adv. Water Resour. 1997, 20, 253–278. [Google Scholar] [CrossRef] [Green Version]
  41. Sun, Y.; Sanaei, P.; Kondic, L.; Cummings, L.J. Modeling and design optimization for pleated membrane filters. Phys. Rev. Fluids 2020, 5, 044306. [Google Scholar] [CrossRef]
  42. Colecchio, I.; Boschan, A.; Otero, A.D.; Noetinger, B. On the multiscale characterization of effective hydraulic conductivity in random heterogeneous media: A historical survey and some new perspectives. Adv. Water Resour. 2020, 140, 103594. [Google Scholar] [CrossRef]
  43. Roa, R.; Zholkovskiy, E.K.; Nägele, G. Ultrafiltration modeling of non-ionic microgels. Soft Matter 2015, 11, 4106–4122. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Probstein, R.F.; Brenner, H. Physicochemical Hydrodynamics: An Introduction; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1994. [Google Scholar]
  45. Van Dyke, M. Perturbation Methods in Fluid Mechanics; The Parabolic Press: Stanford, CA, USA, 1975. [Google Scholar]
  46. Quezada, C.; Estay, H.; Cassano, A.; Troncoso, E.; Ruby-Figueroa, R. Prediction of permeate flux in ultrafiltration processes: A review of modeling approaches. Membranes 2021, 11, 368. [Google Scholar] [CrossRef]
  47. Ghidossi, R.; Veyret, D.; Moulin, P. Computational fluid dynamics applied to membranes: State of the art and opportunities. Chem. Eng. Process. Process. Intensif. 2006, 45, 437–454. [Google Scholar] [CrossRef]
  48. Belfort, G. Membrane filtration with liquids: A global approach with prior successes, new developments and unresolved challenges. Angew. Chem. Int. Ed. 2019, 58, 1892–1902. [Google Scholar] [CrossRef]
  49. Zhou, Z.; Ling, B.; Battiato, I.; Husson, S.M.; Ladner, D.A. Concentration polarization over reverse osmosis membranes with engineered surface features. J. Membr. Sci. 2020, 617, 118199. [Google Scholar] [CrossRef]
  50. Linkhorst, J.; Lölsberg, J.; Thill, S.; Lohaus, J.; Lüken, A.; Nägele, G.; Wessling, M. Templating the morphology of soft microgel assemblies using a nanolithographic 3D-printed membrane. Sci. Rep. 2021, 11, 812. [Google Scholar] [CrossRef]
  51. Linkhorst, J.; Rabe, J.; Hirschwald, L.T.; Kuehne, A.J.C.; Wessling, M. Direct observation of deformation in microgel filtration. Sci. Rep. 2019, 9, 18998. [Google Scholar] [CrossRef] [Green Version]
  52. Bowen, W.R.; Williams, P.M. Quantitative predictive modelling of ultrafiltration processes: Colloidal science approaches. Adv. Colloid Interface Sci. 2007, 134, 3–14. [Google Scholar] [CrossRef] [PubMed]
  53. Heinen, M.; Zanini, F.; Roosen-Runge, F.; Fedunová, D.; Zhang, F.; Hennig, M.; Seydel, T.; Schweins, R.; Sztucki, M.; Antalík, M.; et al. Viscosity and diffusion: Crowding and salt effects in protein solutions. Soft Matter 2012, 8, 1404–1419. [Google Scholar] [CrossRef] [Green Version]
  54. Roa, R.; Menne, D.; Riest, J.; Buzatu, P.; Zholkovskiy, E.K.; Dhont, J.K.G.; Wessling, M.; Nägele, G. Ultrafiltration of charge-stabilized dispersions at low salinity. Soft Matter 2016, 12, 4638–4653. [Google Scholar] [CrossRef] [Green Version]
  55. Brito, M.E.; Denton, A.R.; Nägele, G. Modeling deswelling, thermodynamics, structure, and dynamics in ionic microgel suspensions. J. Chem. Phys. 2019, 151, 224901. [Google Scholar] [CrossRef]
  56. Shen, J.J.S.; Probstein, R.F. On the prediction of limiting flux in laminar ultrafiltration of macromolecular solutions. Ind. Eng. Chem. Fundam. 1977, 16, 459–465. [Google Scholar] [CrossRef]
  57. Probstein, R.F.; Shen, J.S.; Leung, W.F. Ultrafiltration of macromolecular solutions at high polarization in laminar channel flow. Desalination 1978, 24, 1–16. [Google Scholar] [CrossRef]
  58. Romero, C.A.; Davis, R.H. Global model of crossflow microfiltration based on hydrodynamic particle diffusion. Chem. Eng. Sci. 1988, 39, 157–185. [Google Scholar] [CrossRef]
  59. Davis, R.H.; Sherwood, J. A similarity solution for steady-state crossflow microfiltration. Chem. Eng. Sci. 1990, 45, 3203–3209. [Google Scholar] [CrossRef]
  60. Park, G.W. gwpark-git/mBLA_UF: Python code for modified boundary layer solution of concentration-polarization and flow properties in ultrafiltration (v3.1). Zenodo 2021. [Google Scholar] [CrossRef]
Figure 1. Considered membrane geometries: (a) cylindrical membrane of inner radius R (left), and flat sheets membranes of height 2 R and width W R (right). (b) Cross-sectional structure of a cylindrical membrane (left, CM), a membrane consisting of an upper and lower flat sheet part at vertical distance 2 R (middle, FMM), and a flat membrane–flat substrate combination where the bottom substrate sheet is impermeable to particles and solvent (right, FMS). The membrane thickness is h, and L with L R is the axial length of the membrane.
Figure 1. Considered membrane geometries: (a) cylindrical membrane of inner radius R (left), and flat sheets membranes of height 2 R and width W R (right). (b) Cross-sectional structure of a cylindrical membrane (left, CM), a membrane consisting of an upper and lower flat sheet part at vertical distance 2 R (middle, FMM), and a flat membrane–flat substrate combination where the bottom substrate sheet is impermeable to particles and solvent (right, FMS). The membrane thickness is h, and L with L R is the axial length of the membrane.
Membranes 11 00960 g001
Figure 2. (a) Single-particle shear-Pèclet number, P e a a 3 , as a function of particle radius a and for a shear rate γ ˙ * = 270 / s typical of UF. (b) Carnahan–Starling-based osmotic pressure at freezing, Π ( ϕ f , a ) , and thermal pressure, k B T / V a , as functions of particle radius a. Open circles indicate the radii a = 3.13 nm and a = 10 nm used in this work.
Figure 2. (a) Single-particle shear-Pèclet number, P e a a 3 , as a function of particle radius a and for a shear rate γ ˙ * = 270 / s typical of UF. (b) Carnahan–Starling-based osmotic pressure at freezing, Π ( ϕ f , a ) , and thermal pressure, k B T / V a , as functions of particle radius a. Open circles indicate the radii a = 3.13 nm and a = 10 nm used in this work.
Membranes 11 00960 g002
Figure 3. Reduced osmotic pressure, Π / ( n k B T ) , as a function of particle volume fraction ϕ . The solid line is the Carnahan–Starling prediction for hard spheres, compared to experimental data (open symbols) for bovine serum albumin solutions at different pH values as indicated. Experimental data are reproduced from [30]. Dashed lines are empirical fits to the data.
Figure 3. Reduced osmotic pressure, Π / ( n k B T ) , as a function of particle volume fraction ϕ . The solid line is the Carnahan–Starling prediction for hard spheres, compared to experimental data (open symbols) for bovine serum albumin solutions at different pH values as indicated. Experimental data are reproduced from [30]. Dashed lines are empirical fits to the data.
Membranes 11 00960 g003
Figure 4. (a) Concentration dependence of the short-time sedimentation coefficient K s e d ( ϕ ) . Solid line is the prediction by Equation (4), and open symbols are dynamic simulation data by three groups as indicated [32,33,34]. (b) Short-time collective diffusion coefficient, D ( ϕ ) , according to Equation (3) (solid line), given in units of D 0 . (c) Dispersion shear viscosity, η ( ϕ ) , in units of the solvent viscosity η s . Solid line is the prediction by Equations (5)–(7). Open symbols are dynamic simulation and filled symbols experimental viscosity data reproduced from Refs. [35,36] and Refs. [33,37], respectively.
Figure 4. (a) Concentration dependence of the short-time sedimentation coefficient K s e d ( ϕ ) . Solid line is the prediction by Equation (4), and open symbols are dynamic simulation data by three groups as indicated [32,33,34]. (b) Short-time collective diffusion coefficient, D ( ϕ ) , according to Equation (3) (solid line), given in units of D 0 . (c) Dispersion shear viscosity, η ( ϕ ) , in units of the solvent viscosity η s . Solid line is the prediction by Equations (5)–(7). Open symbols are dynamic simulation and filled symbols experimental viscosity data reproduced from Refs. [35,36] and Refs. [33,37], respectively.
Membranes 11 00960 g004
Figure 5. (a) mBLA prediction for the longitudinal velocity factor U ( y , z , [ ϕ w ] ) , using ϕ w = 0.4 , a = 10 nm, ϵ δ 1.28 × 10 2 , and v w = 3.35 × 10 6 m/s. The inset shows the stretched U / ϵ δ as a function of the stretched distance from the bottom wall, ( y + R ) / δ C P . While the bottom wall is a permeable membrane for FMM, it is impermeable for FMS. (b) Transversal velocity factor V ( y ) = V o u t ( y ) in Equation (31) for CM, FMM, and FMS, respectively.
Figure 5. (a) mBLA prediction for the longitudinal velocity factor U ( y , z , [ ϕ w ] ) , using ϕ w = 0.4 , a = 10 nm, ϵ δ 1.28 × 10 2 , and v w = 3.35 × 10 6 m/s. The inset shows the stretched U / ϵ δ as a function of the stretched distance from the bottom wall, ( y + R ) / δ C P . While the bottom wall is a permeable membrane for FMM, it is impermeable for FMS. (b) Transversal velocity factor V ( y ) = V o u t ( y ) in Equation (31) for CM, FMM, and FMS, respectively.
Membranes 11 00960 g005
Figure 6. CM, FMM, and FMS concentration profiles at the membrane wall, ϕ w ( z ) , for a dispersion of Brownian hard spheres of radius a = 3.13 nm (a) and a = 10 nm (b), respectively. The insets show the according axial velocity profiles, u ( y = 0 , z ) , at the channel center-line y = 0 , in units of the CM inlet velocity u C M 0 = u ( 0 , 0 ) . Black curves are for the same mean inlet velocity u ¯ 0 . Dashed (red) curves represent a second FMM system (referred to as FMM2), having a flow shear rate at the membrane inlet equal to that for the CM geometry.
Figure 6. CM, FMM, and FMS concentration profiles at the membrane wall, ϕ w ( z ) , for a dispersion of Brownian hard spheres of radius a = 3.13 nm (a) and a = 10 nm (b), respectively. The insets show the according axial velocity profiles, u ( y = 0 , z ) , at the channel center-line y = 0 , in units of the CM inlet velocity u C M 0 = u ( 0 , 0 ) . Black curves are for the same mean inlet velocity u ¯ 0 . Dashed (red) curves represent a second FMM system (referred to as FMM2), having a flow shear rate at the membrane inlet equal to that for the CM geometry.
Membranes 11 00960 g006
Figure 7. Excess and bulk axial flux parts, J z e x ( y , z ) (solid lines) and J z b ( y , z ) (dashed lines), plotted as functions of the reduced distance, 1 y / R , from the membrane wall, for systems CM (left) and FMM (right) using a = 3.13 nm. The fluxes are given in units of J z ¯ 0 = ϕ b u ¯ 0 . The arrows mark increasing values of z, with z / L = 0.2 , 0.4 , 0.6 , 0.8 , 1.0 . The dotted vertical lines mark the distance from the membrane wall equal to δ C P . System parameters as in Figure 6a.
Figure 7. Excess and bulk axial flux parts, J z e x ( y , z ) (solid lines) and J z b ( y , z ) (dashed lines), plotted as functions of the reduced distance, 1 y / R , from the membrane wall, for systems CM (left) and FMM (right) using a = 3.13 nm. The fluxes are given in units of J z ¯ 0 = ϕ b u ¯ 0 . The arrows mark increasing values of z, with z / L = 0.2 , 0.4 , 0.6 , 0.8 , 1.0 . The dotted vertical lines mark the distance from the membrane wall equal to δ C P . System parameters as in Figure 6a.
Membranes 11 00960 g007
Figure 8. Cross-section averaged axial excess particle-flux, J z ¯ e x ( z ) , in units of J z ¯ 0 , as function of reduced axial distance, z / L , from the inlet. Two dispersions with particle radius a = 3.13 nm (a) and a = 10 nm (b) are considered for input parameters as in Figure 6a,b, respectively.
Figure 8. Cross-section averaged axial excess particle-flux, J z ¯ e x ( z ) , in units of J z ¯ 0 , as function of reduced axial distance, z / L , from the inlet. Two dispersions with particle radius a = 3.13 nm (a) and a = 10 nm (b) are considered for input parameters as in Figure 6a,b, respectively.
Membranes 11 00960 g008
Figure 9. (a-1,b-1): TMP dependence of average permeate flux v w , solvent recovery indicator β , and concentration factor α for particle radius a = 3.13 nm (panels a-1a-3) and a = 10 nm (panels b-1b-3), respectively. Symbols are mBLA results, and solid lines are pure solvent predictions. Input parameters except for TMP are as in Figure 6.
Figure 9. (a-1,b-1): TMP dependence of average permeate flux v w , solvent recovery indicator β , and concentration factor α for particle radius a = 3.13 nm (panels a-1a-3) and a = 10 nm (panels b-1b-3), respectively. Symbols are mBLA results, and solid lines are pure solvent predictions. Input parameters except for TMP are as in Figure 6.
Membranes 11 00960 g009
Figure 10. Reduced average permeate flux, v w ( R / D 0 ) , (left) and concentration factor, α , (right) versus feed concentration, ϕ b , for fixed P e R 78 . Input parameters are as in Figure 6, except for ϕ b , which is varied. Open symbols are for a = 3.13 nm and filled symbols for a = 10 nm. The inset in the left panel shows the pressure ratio Π / TMP . The plateau values, α 0 , of α are marked in the right panel by the horizontal dashed line segments.
Figure 10. Reduced average permeate flux, v w ( R / D 0 ) , (left) and concentration factor, α , (right) versus feed concentration, ϕ b , for fixed P e R 78 . Input parameters are as in Figure 6, except for ϕ b , which is varied. Open symbols are for a = 3.13 nm and filled symbols for a = 10 nm. The inset in the left panel shows the pressure ratio Π / TMP . The plateau values, α 0 , of α are marked in the right panel by the horizontal dashed line segments.
Membranes 11 00960 g010
Figure 11. Solvent recovery indicator, β , as function of the cross-section averaged inlet velocity u ¯ 0 in units of L D 0 / R 2 . Here, u ¯ 0 is varied for a fixed TMP = 16 kPa for a = 3.13 nm, and 5 kPa for a = 10 nm. The feed concentration is ϕ b = 0.001 . Remaining input parameters as in Figure 6a,b. Open (filled) symbols are mBLA results for a = 3.13 nm ( a = 10 nm). Solid, dashed, and dotted lines represent the pure solvent values, β 0 , for the respective geometries as indicated.
Figure 11. Solvent recovery indicator, β , as function of the cross-section averaged inlet velocity u ¯ 0 in units of L D 0 / R 2 . Here, u ¯ 0 is varied for a fixed TMP = 16 kPa for a = 3.13 nm, and 5 kPa for a = 10 nm. The feed concentration is ϕ b = 0.001 . Remaining input parameters as in Figure 6a,b. Open (filled) symbols are mBLA results for a = 3.13 nm ( a = 10 nm). Solid, dashed, and dotted lines represent the pure solvent values, β 0 , for the respective geometries as indicated.
Membranes 11 00960 g011
Figure 12. Solvent recovery indicator, β , as a function of P e R , for ( β 0 / P e R ) G 1 , G 2 8.08 × 10 3 , and ( γ 0 / P e R ) G 1 1.55 × 10 4 (list G1 in Table 3: open symbols) and ( γ 0 / P e R ) G 2 3.10 × 10 4 (list G2: filled symbols), respectively. Identical membrane properties are used for CM, FMM, and FMS. The employed operating parameters are listed in rows G1 and G2 of Table 3.
Figure 12. Solvent recovery indicator, β , as a function of P e R , for ( β 0 / P e R ) G 1 , G 2 8.08 × 10 3 , and ( γ 0 / P e R ) G 1 1.55 × 10 4 (list G1 in Table 3: open symbols) and ( γ 0 / P e R ) G 2 3.10 × 10 4 (list G2: filled symbols), respectively. Identical membrane properties are used for CM, FMM, and FMS. The employed operating parameters are listed in rows G1 and G2 of Table 3.
Membranes 11 00960 g012
Figure 13. Solvent recovery indicator, β , as a function of P e R for three different variable sets ( β 0 / P e R , γ 0 / P e R ) given in the figure for the respective CM, FMM, and FMS geometries. Symbols are mBLA results for β based on the according input parameter lists in Table 3. See text for details.
Figure 13. Solvent recovery indicator, β , as a function of P e R for three different variable sets ( β 0 / P e R , γ 0 / P e R ) given in the figure for the respective CM, FMM, and FMS geometries. Symbols are mBLA results for β based on the according input parameter lists in Table 3. See text for details.
Membranes 11 00960 g013
Table 1. Summary of geometry-dependent quantities characterizing the outer solution. Notice that λ 2 = R H M / ( L A U ¯ o u t ) . The effective permeability parameter, K, is given in units of K * = η s L p L 2 / R 3 so that K / K * = λ 1 λ 2 ( R / R H ) 3 .
Table 1. Summary of geometry-dependent quantities characterizing the outer solution. Notice that λ 2 = R H M / ( L A U ¯ o u t ) . The effective permeability parameter, K, is given in units of K * = η s L p L 2 / R 3 so that K / K * = λ 1 λ 2 ( R / R H ) 3 .
Membrane GeometryTransversal Velocity
Boundary Condition
λ 1 λ 2 R H U ¯ out M A K K * H V out ( y )
CM v ( R , z ) = v w ( z ) 12 R 2 1 2 2 L R 4 R ln 1 + h R 2 y R y R 3
FMM v ( R , z ) = v w ( z ) v ( R , z ) = v w ( z ) 2 3 2 R 2 3 L R 3 h 1 2 3 y R y R 3
FMS v ( R , z ) = v w ( z ) v ( R , z ) = 0 2 3 4 R 2 3 L 2 R 3 2 h 1 4 3 y R y R 3 + 2
Table 2. Summary of conditions for the validity of the mBLA method. The Reynolds number associated with the channel flow was R e = 4 R H u ¯ 0 ρ / η , where ρ is the dispersion mass density, η the effective dispersion viscosity, u ¯ 0 the cross-section averaged inlet velocity, and R H the hydraulic channel radius. The solvent recovery indicator, β , is defined in Equation (21).
Table 2. Summary of conditions for the validity of the mBLA method. The Reynolds number associated with the channel flow was R e = 4 R H u ¯ 0 ρ / η , where ρ is the dispersion mass density, η the effective dispersion viscosity, u ¯ 0 the cross-section averaged inlet velocity, and R H the hydraulic channel radius. The solvent recovery indicator, β , is defined in Equation (21).
ConditionsRemarks
P e a 0.1 Strong Brownian motion
R / L 1 Small aspect ratio. Note that ϵ δ R / L
R e 2000 Condition for laminar (non-turbulent) flow
R e R / L 1 No inertial flow effects on length scale L
β O ϵ δ Condition for significant permeability effects
ϕ b 1 Small feed concentration (i.e., ϕ b P e R < 1 )
Table 3. Input parameters used for the mBLA results for the solvent recovery indicator β depicted in Figure 12 and Figure 13, respectively, as indicated. Fixed parameters are particle radius a = 3.13 nm and mean Darcy permeability κ = 1.36 × 10 16 m 2 . The employed values for the reference system cross-section averaged inlet velocity and hydraulic permeability are u ¯ R E F 0 = 3.40 × 10 2 and L p , R E F = 6.7 × 10 10 m/(Pa s), respectively.
Table 3. Input parameters used for the mBLA results for the solvent recovery indicator β depicted in Figure 12 and Figure 13, respectively, as indicated. Fixed parameters are particle radius a = 3.13 nm and mean Darcy permeability κ = 1.36 × 10 16 m 2 . The employed values for the reference system cross-section averaged inlet velocity and hydraulic permeability are u ¯ R E F 0 = 3.40 × 10 2 and L p , R E F = 6.7 × 10 10 m/(Pa s), respectively.
Data
(Figure 12)
R(mm)L(m) h / R L p / L p , REF
CM, FMM/FMS
ϕ b / 10 3
CM, FMM, FMS
u ¯ 0 / u ¯ REF 0
CM, FMM, FMS
( β 0 / Pe R ) / 10 3 ( γ 0 / Pe R ) / 10 4
G10.50.5-1.00, 1.001.00, 0.375, 0.1861.00, 0.50, 0.258.081.55
G20.50.5-1.00, 1.002.00, 0.750, 0.3751.00, 0.50, 0.258.083.10
(Figure 13) R (mm) L (m) h / R L p / L p , REF ϕ b / 10 3 u ¯ 0 / u ¯ REF 0 ( β 0 / Pe R ) / 10 3 ( γ 0 / Pe R ) / 10 4
CM-10.50.5-1.001.001.008.081.55
CM-20.50.50.51.001.001.008.081.55
CM-30.50.510.591.691.008.081.55
CM-40.250.50.52.001.004.008.081.55
FMM-10.50.5-1.001.001.004.042.07
FMM-20.50.50.50.811.241.004.042.08
FMM-30.50.510.412.441.004.042.07
FMM-40.250.50.51.621.234.004.042.07
FMS-10.50.5-1.001.001.002.022.07
FMS-20.50.50.50.811.241.002.022.08
FMS-30.50.510.412.441.002.022.07
FMS-40.250.50.51.621.234.002.022.07
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Park, G.W.; Nägele, G. Geometrical Influence on Particle Transport in Cross-Flow Ultrafiltration: Cylindrical and Flat Sheet Membranes. Membranes 2021, 11, 960. https://doi.org/10.3390/membranes11120960

AMA Style

Park GW, Nägele G. Geometrical Influence on Particle Transport in Cross-Flow Ultrafiltration: Cylindrical and Flat Sheet Membranes. Membranes. 2021; 11(12):960. https://doi.org/10.3390/membranes11120960

Chicago/Turabian Style

Park, Gun Woo, and Gerhard Nägele. 2021. "Geometrical Influence on Particle Transport in Cross-Flow Ultrafiltration: Cylindrical and Flat Sheet Membranes" Membranes 11, no. 12: 960. https://doi.org/10.3390/membranes11120960

APA Style

Park, G. W., & Nägele, G. (2021). Geometrical Influence on Particle Transport in Cross-Flow Ultrafiltration: Cylindrical and Flat Sheet Membranes. Membranes, 11(12), 960. https://doi.org/10.3390/membranes11120960

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop