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

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 . 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 and width W, where and . 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 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.

Figure 1

Considered membrane geometries: (a) cylindrical membrane of inner radius R (left), and flat sheets membranes of height and width (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 (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 is the axial length of the membrane.

The CM membrane is distinguished from the FMM membrane by a larger cross-sectional circumference-to-area ratio of , while this ratio was equal to in the FMM case. For otherwise equal membrane properties, i.e., equal hydraulic permeability , 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 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, , obeys [26,27] 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, is the single-particle Stokes–Einstein diffusion coefficient for stick hydrodynamic boundary conditions [26], and is the viscosity of the Newtonian solvent. For small , 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.

Figure 2

(a) Single-particle shear-Pèclet number, , as a function of particle radius a and for a shear rate typical of UF. (b) Carnahan–Starling-based osmotic pressure at freezing, , and thermal pressure, , as functions of particle radius a. Open circles indicate the radii nm and nm used in this work.

The equilibrium osmotic pressure, , of Brownian hard spheres can be described by the Carnahan–Starling expression [28,29] where is the particle volume, and is the particle volume fraction. This expression holds to high accuracy up to the freezing transition volume fraction . In Figure 2b, the maximal fluid-state osmotic pressure reached at freezing, , is plotted as a function of the particle radius a in conjunction with the thermal pressure , where is the Boltzmann constant, and T is the dispersion temperature. Notice that the 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 in a M saline solution. We determined the volume fractions in Figure 3 from the protein concentration values and the equivalent spherical radius value nm given in [30]. The pressure data at are located below the Carnahan–Starling curve, suggesting that there is an attractive interaction contribution.

Figure 3

Reduced osmotic pressure, , 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.

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] where is the short-time sedimentation coefficient, is the single-particle diffusion coefficient, and is the osmotic compressibility factor. For the sedimentation coefficient of Brownian hard spheres, we use the analytic expression [31]

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, , the accurate Carnahan–Starling-type expression 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, , should not be confused with the self-diffusion coefficient, . In fact, the latter coefficient decreases with increasing concentration below its infinite dilution value, , both for attractive and repulsive interactions.

Figure 4

(a) Concentration dependence of the short-time sedimentation coefficient . 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, , according to Equation (3) (solid line), given in units of . (c) Dispersion shear viscosity, , in units of the solvent viscosity . 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.

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 , the two viscosity contributions are approximated, to good accuracy up to the freezing volume fraction, by [31] with the Saito-type function . Here, is a good analytic approximation for the short-time self-diffusion coefficient of Brownian hard spheres. Up to a factor of six, 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, , is accordingly of the form where is the unit vector in axial direction. For the CM channel, is the radial unit vector pointing from the inner axis to the cylindrical membrane wall, whereas for the FMM and FMS channels, is the unit vector in the transversal direction, perpendicular to the parallel sheets (c.f. Figure 1). Furthermore, is the transversal, and 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 of particles at position and time t, and the dispersion-averaged velocity field .

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], in conjunction with the incompressibility constraint, , for the dispersion velocity field . The dispersion-averaged pressure field, , is the sum of the equilibrium particle osmotic pressure, , and a fluid-phase pressure contribution, , adjusting itself such that incompressibility is fulfilled. There is a contribution to the Stokes equation proportional to /, operative in the inhomogeneous CP layer region.

The dispersion-averaged, stationary local particle flux, , with has an advection contribution, , and a diffusion contribution, , respectively, with the latter quantified by the gradient diffusion coefficient introduced in Equation (3). Substituting this flux into the macroscopic continuity equation expressing mass conservation leads to the steady-state advection–diffusion equation, where . Since we assume fully particle-retentive membranes, it holds that 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] for the transversal flow velocity (permeate flux), , at the inner membrane wall(s). The permeate flux is here defined with a positive sign. In FMS geometry, is obviously zero at the impermeable lower substrate wall. Here, is the local transmembrane pressure at axial distance z from the inlet, and and are the lumen-side pressure and particle volume concentration at the membrane wall(s), respectively. The pressure at the permeate side, , is taken as constant.

The hydraulic permeability, , of a clean membrane is given by [38,39] 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 , and Here, is the feed concentration of particles at the inlet, and is the pressure given at the outlet. The velocity field at the inlet port (where ) is taken here as a fully developed parabolic flow in axial direction, so that . Furthermore, is the axial velocity at the inlet center.

The present UF model is specified by the input (operating) parameters . The pressure at the outlet is fixed to 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 , , and atm, and instead of specifying the permeate pressure, we use alternatively as a fourth input parameter the channel-length-averaged linearized transmembrane presssure, , 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., associated with a Hagen–Poiseuille-type quadratic velocity field inside the CM, FMM, and FMS channels occurring for and . The brackets denote the channel-length average,

The linearized TMP and the actual TMP, , are thus expressed as respectively. Here, is the, in general, non-linear dispersion pressure profile in UF. Moreover, it is , with the hydraulic radius and the dimensionless geometry coefficient for , respectively (see Table 1). The reason for using the linearized TMP as input parameter is that for values of 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 .

Table 1

Summary of geometry-dependent quantities characterizing the outer solution. Notice that . The effective permeability parameter, K, is given in units of so that .

Membrane Geometry	Transversal Velocity Boundary Condition	λ_1	λ_2	R_H	U¯^out	MA	KK^*	H	V^out(y)
CM	v(R,z)=v_w(z)	1	2	R2	12	2LR	4	Rln1+hR	2yR−yR^3
FMM	v(R,z)=v_w(z)v(−R,z)=−v_w(z)	2	32	R	23	LR	3	h	123yR−yR^3
FMS	v(R,z)=v_w(z)v(−R,z)=0	2	34	R	23	L2R	32	h	143yR−yR^3+2

A useful process indicator characterizing the filtration efficiency is the solvent recovery indicator, , given by [2,43] Here, is the cross-sectional average of the inlet velocity . The overline indicates the cross-sectional average, with the constant cross-sectional area equal to for CM, and both for FMM and FMS. Furthermore, is the dispersion volume flow rate through the inlet cross-section, and is the permeate volume flow rate through the inner membrane area M. This area is equal to for , respectively. A larger value for reflects a larger concentration of particles at the outlet port. Notice that volume conservation implies that , where 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 , across which is steeply decreasing, from its maximal value, , attained at the wall towards its minimal bulk value . On introducing the smallness parameter with , 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, , i.e., [15,25,44] The transversal Pèclet number is the ratio of the diffusion time of particles, , and the transversal flow advection time, , across the transversal distance R. It is hereby assumed that the feed concentration at the inlet is small, i.e., . For a significantly developed CP layer, is of the order . 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] 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 and with longitudinal velocity factor and transversal velocity factor

To zeroth order in , the outer velocity field is factorized in y and z. The permeate flux, depends on the pressure field, , through the Darcy–Starling expression, while the axial velocity at the center-line, , is determined by the pressure according to .

The pressure, in turn, is obtained as where and We have introduced here the dimensionless effective permeability parameter, K, given by 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 and , with for , on the (curvature-corrected) membrane thickness given by for CM, and for FMM and FMS, respectively, and on the cross-sectional average, , of equal to for CM and 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 . The leading-order continuity, effective Stokes, and advection–diffusion equations determining the inner solutions are obtained as [15] respectively. The first (continuity) and second (Stokes) equation state that the transversal velocity, v, and the shear stress, , are independent of y inside the CP layer, to leading order in . Thus, 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 where the functional dependence on the (up to this point) unknown concentration profile, , is indicated. We have introduced the dimensionless function appearing in the inner solutions for and u. The transversal variable y is restricted to the interval for the FMS geometry, which lacks reflection symmetric with respect to the plane . For CM and FMM, y is further restricted to , since, for FMM, it holds that 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 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 The fields and 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 with . For constant viscosity, it reduces to , where is the parabolic profile in Equation (30).

The remaining task is to determine the particle concentration profile, , at the membrane wall. To this end, we employ, as a global condition, the cross-sectional particle-flux conservation law, where 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 and of assuming a fully particle retentive membrane. We ignore here a small longitudinal diffusion flux contribution of , with the implication that . Notice that is determined, according to the inlet boundary conditions, by and . We deviate here from our earlier mBLA work in [15] where the inlet pressure was prescribed as an inlet boundary value instead of .

Using particle-flux conservation combined with a fixed-point iteration method adapted to the flat sheets FMM and FMS geometries, is numerically determined. The concentration and flow fields follow then from substituting 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, 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 . Differences from the parabolic profile are visible in the CP boundary layer region, provided they are enhanced by dividing through . Figure 5a shows a simplified mBLA result for , obtained by neglecting for simplicity the z-dependence of the matched solutions in Equations (43)–(46), and using constant , , nm, , and mm. For the suspension transport properties 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 except for the boundary layer region magnified in the inset. In the inset, is plotted as a function of the stretched distance, , from the bottom wall (in case of FMM and FMS). Regarding the FMS geometry, a slight deviation from symmetry with respect to the 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.

Figure 5

(a) mBLA prediction for the longitudinal velocity factor , using , nm, , and m/s. The inset shows the stretched as a function of the stretched distance from the bottom wall, . While the bottom wall is a permeable membrane for FMM, it is impermeable for FMS. (b) Transversal velocity factor in Equation (31) for CM, FMM, and FMS, respectively.

For comparison, in Figure 5b, the velocity factor 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, for CM has its peak value attained away from the membrane wall at radial distance from the center-line. While in CM and FMM due to symmetry, the curve of in FMS is strongly asymmetric with , 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.

Table 2

Summary of conditions for the validity of the mBLA method. The Reynolds number associated with the channel flow was , where is the dispersion mass density, the effective dispersion viscosity, the cross-section averaged inlet velocity, and the hydraulic channel radius. The solvent recovery indicator, , is defined in Equation (21).

Conditions	Remarks
Pe_a≤0.1	Strong Brownian motion
R/L≪1	Small aspect ratio. Note that ϵ_δ≥R/L
Re≲2000	Condition for laminar (non-turbulent) flow
ReR/L≪1	No inertial flow effects on length scale L
β≥Oϵ_δ	Condition for significant permeability effects
ϕ_b≪1	Small feed concentration (i.e., ϕ_bPe_R<1)

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 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, , axial velocity at the inlet center, , and linearized TMP, . The pressure at the outlet is fixed to atm. These parameters are complemented by the selected membrane and dispersion parameters including R and L, membrane hydraulic permeability, , solvent viscosity, , 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 , mm, m, and m/(Pa s). The selected value for is in the range of typical UF permeabilities [2], and nm. Reference values for depend on the membrane geometry and are selected as m/s for , respectively. In all three geometries, the same value, m/s, is obtained for the cross-section averaged inlet velocity, where and . 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, 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 , the geometry-dependent reference values for the effective permeability parameter K defined in Equation (35) are obtained as for , respectively. Since , the actual TMP, , is practically equal to the linearized . 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, , and for the longitudinal particle flux . For the latter quantity, both its bulk and excess parts are discussed. The TMP is not varied here, but two different particle radii 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, , 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, , at the membrane wall for particle radii (a) nm and (b) nm, respectively, for fixed transversal Pèclet number , and for the reference values of , R, and L. Consider first Figure 6a for the concentration profile, where the reference value for the average inlet velocity, 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 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 , can be qualitatively understood from inspecting (in the inset) the axial velocity, , at the axial center-line . Note that for the CM geometry is larger than for FMM and FMS, by a factor of . At axial distances similar to those for the profile, crossovers are observed between the three curves for . Since , a stronger decay of reflects a larger transmembrane flux, , and thus an enhanced transport of particles towards the membrane wall.

Figure 6

CM, FMM, and FMS concentration profiles at the membrane wall, , for a dispersion of Brownian hard spheres of radius nm (a) and nm (b), respectively. The insets show the according axial velocity profiles, , at the channel center-line , in units of the CM inlet velocity . Black curves are for the same mean inlet velocity . 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.

The usage of the same value for in the three geometries implies that the shear rate at the inlet of the membrane wall, (c.f. Equation (1)), is larger for CM than for FMM and FMS. This explains the initially smaller increase in 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 at the inlet is the same for CM and FMM2, the values for FMM2 sufficiently distant from the inlet are distinctly smaller than in the CM case. This is reflected in the axial velocity profile 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 nm, smaller TMP kPa, and cross-section averaged inlet velocity 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 at the channel outlet.

The ratio, , of length-averaged transmembrane osmotic pressure and TMP equals and in Figure 6a,b, respectively. It reflects the strong influence of the CP layer on the permeate flux , 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 can be non-monotonic in z, with maximum concentration attained for . As discussed in [15], in the context of the CM geometry, a criterion for a strictly monotonic increase in is that .

We proceed by discussing the axial particle flux, , which according to consists of an excess flux contribution, , that vanishes in the bulk region and a bulk flux contribution, . Right at the membrane wall, due to the employed stick boundary condition. Without a significant CP layer, is practically zero. The results by the mBLA method for the excess and bulk parts of are depicted in Figure 7 for CM and FMM, respectively, using nm and input parameters as in Figure 6a. The fluxes are plotted semi-logarithmically as functions of the (relative) distance, , 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, , 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 ( towards its maximal value at . The characteristic CP layer thickness, , 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 is located at a distance from the wall larger than , 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 at the channel center-line .

Figure 7

Excess and bulk axial flux parts, (solid lines) and (dashed lines), plotted as functions of the reduced distance, , from the membrane wall, for systems CM (left) and FMM (right) using nm. The fluxes are given in units of . The arrows mark increasing values of z, with , . The dotted vertical lines mark the distance from the membrane wall equal to . System parameters as in Figure 6a.

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 , reflected in the likewise different behavior of 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, , is shown for nm and nm, respectively. Like , 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 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 Recall next from Table 1 that the ratio, , 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 , which arises since the value for in FMS is half that of FMM, while 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.

Figure 8

Cross-section averaged axial excess particle-flux, , in units of , as function of reduced axial distance, , from the inlet. Two dispersions with particle radius nm (a) and nm (b) are considered for input parameters as in Figure 6a,b, respectively.

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

Having discussed the local profiles and , 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, , and the solvent recovery indicator 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, , defined as the ratio of cross-section averaged concentration at the outlet and the (cross-section averaged) feed concentration . Owing to particle flux and dispersion volume conservation valid for an ideally particle retentive membrane, the concentration factor is related to simply by [43] One generally observes that 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 , , 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 are nearly coincidental, even for the largest TMP considered. The straight solid lines depict the pure-solvent form , for the same hydraulic permeability used in both figure panels. The influence of the CP layer (osmotic pressure) buildup on the permeate flux is noticeable at larger TMP where 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 kPa.

Figure 9

(a-1,b-1): TMP dependence of average permeate flux , solvent recovery indicator , and concentration factor for particle radius nm (panels a-1–a-3) and nm (panels b-1–b-3), respectively. Symbols are mBLA results, and solid lines are pure solvent predictions. Input parameters except for TMP are as in Figure 6.

The straight lines in panels (a-2) and (b-2) of Figure 9 are the pure-solvent prediction, , equal to at zero feed concentration . Explicitly, where M is the lumen-side membrane surface area, and A is the area of the inner cross-section. As noticed from Equation (51), has a geometry-dependent factor, , of values listed in Table 1, and it includes the cross-section averaged inlet velocity , which has different values in panels (a) and (b), respectively. At small TMP, the mBLA curves for overlap with the pure-solvent lines of , since the osmotic pressure contribution is here negligible. For comparison, panels (a-3) and (b-3) show the associated concentration factor . Solid lines represent , which quantifies the initial increase in with increasing TMP. Notice that, according to Equation (51), both and are solely dependent on given input parameters. According to Table 1, it holds that = 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 , and it is plotted versus multiplied by , i.e., the selected value for the transversal Pèclet number. The multiplication of by is used here to highlight that , an inequality required for the validity of the mBLA method (c.f. Table 2). Except for , all input parameters are as in Figure 6a,b, respectively. Regarding Figure 10 (left), one notices that for . According to the inset, the initial plateau value of the reduced flux, equal to , is mirrored by the plateau value . With further increasing , the accordingly increasing causes a decline in the reduced permeate flux linear in the logarithm of , 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 is close to one. The mBLA results for as a function of are displayed in Figure 10 (right). The concentration factor decreases with increasing , 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.

Figure 10

Reduced average permeate flux, , (left) and concentration factor, , (right) versus feed concentration, , for fixed . Input parameters are as in Figure 6, except for , which is varied. Open symbols are for nm and filled symbols for nm. The inset in the left panel shows the pressure ratio . The plateau values, , of are marked in the right panel by the horizontal dashed line segments.

The dependence of on the cross-section averaged axial inlet velocity, , 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 . With increasingly large inlet velocity, the mBLA data for (symbols) converge towards the pure solvent prediction (lines). This convergence arises since the influence of the CP layer ceases with increasing axial flow velocity.

Figure 11

Solvent recovery indicator, , as function of the cross-section averaged inlet velocity in units of . Here, is varied for a fixed kPa for nm, and 5 kPa for nm. The feed concentration is . Remaining input parameters as in Figure 6a,b. Open (filled) symbols are mBLA results for nm ( nm). Solid, dashed, and dotted lines represent the pure solvent values, , for the respective geometries as indicated.

5.3. Universal Behavior of Global UF Indicators

In Section 5.2, we analyzed the global UF indicators , , and , as functions of TMP, , a, and , 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 , 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, , 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 Here, is the pure solvent recovery indicator introduced in Equation (51) where its linear dependence on TMP was emphasized. While is independent of dispersion properties, these are accounted for by and . The third variable, , includes a geometry-dependent factor for , respectively. Moreover, it depends on the feed concentration , particle radius a, hydraulic membrane permeability , solvent viscosity , and the inner channel half-height R. The appearance and the form of are motivated by a simplified mBLA calculation of presented in Appendix A.

On expressing in terms of the three variables, we find that 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 , for values of the duplet given by the two sets (open symbols) and (filled symbols), respectively. Within each set, the mBLA data for coincide practically for all and for the three geometries CM, FMM, and FMS. The two curves and traced out by the mBLA data for merge at small with the straight dotted line for determined by Equation (52). The two sets G1 and G2 for share the same value for . 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) and , a larger value of causes a smaller values of due to the lowered permeate flux. According to Equation (53), increases for increasing and , and for decreasing particle radius a.

Figure 12

Solvent recovery indicator, , as a function of , for , and (list G1 in Table 3: open symbols) and (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.

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 nm and mean Darcy permeability . The employed values for the reference system cross-section averaged inlet velocity and hydraulic permeability are and 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¯REF0 CM, FMM, FMS	(β^0/Pe_R)/10^−3	(γ^0/Pe_R)/10^−4
G1	0.5	0.5	-	1.00, 1.00	1.00, 0.375, 0.186	1.00, 0.50, 0.25	8.08	1.55
G2	0.5	0.5	-	1.00, 1.00	2.00, 0.750, 0.375	1.00, 0.50, 0.25	8.08	3.10
(Figure 13)	R (mm)	L (m)	h/R	L_p/L_p,REF	ϕ_b/10^−3	u¯^0/u¯REF0	(β^0/Pe_R)/10^−3	(γ^0/Pe_R)/10^−4
CM-1	0.5	0.5	-	1.00	1.00	1.00	8.08	1.55
CM-2	0.5	0.5	0.5	1.00	1.00	1.00	8.08	1.55
CM-3	0.5	0.5	1	0.59	1.69	1.00	8.08	1.55
CM-4	0.25	0.5	0.5	2.00	1.00	4.00	8.08	1.55
FMM-1	0.5	0.5	-	1.00	1.00	1.00	4.04	2.07
FMM-2	0.5	0.5	0.5	0.81	1.24	1.00	4.04	2.08
FMM-3	0.5	0.5	1	0.41	2.44	1.00	4.04	2.07
FMM-4	0.25	0.5	0.5	1.62	1.23	4.00	4.04	2.07
FMS-1	0.5	0.5	-	1.00	1.00	1.00	2.02	2.07
FMS-2	0.5	0.5	0.5	0.81	1.24	1.00	2.02	2.08
FMS-3	0.5	0.5	1	0.41	2.44	1.00	2.02	2.07
FMS-4	0.25	0.5	0.5	1.62	1.23	4.00	2.02	2.07

In rows G1 and G2 of Table 3, the membrane parameters R, L, and 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 , for three sets of the duplet given by (circles), (squares), and (triangles). These three sets are not taken out of the blue but are precisely the values for and of the reference system introduced in Section 5 for CM, FMM, and FMS, respectively. For each geometry with its according set for , 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 , 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 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 for . It holds that (see also Figure 9(a-2)). The mBLA data for in Figure 13 nicely reconfirm the validity of Equation (54).

Figure 13

Solvent recovery indicator, , as a function of for three different variable sets 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.

The universal dependence of on the three variables , , and in Equation (54) is inherited by the concentration factor, , and by the reduced permeate flux, , since these are related to by

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 Here, g is a complicated function given in Equation (A6) of the Appendix A, which depends on 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 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 , , and are the reduced set of independent variables and that it provides, in addition, the expression for 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 . 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: , , and . 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.