Filtration with Multiple Species of Particles

Membrane filtration of feed containing multiple species of particles is a common process in the industrial setting. In this work, we propose a model for filtration of a suspension containing multiple particle species (concrete examples of our model are shown in two and three species), each with different affinities for the material of the porous filter membrane. Using the pore shape within the membrane as a design objective, we formulate a number of optimization problems pertaining to effective separation of desired and undesired particles in the special case of two-particle species and we present results showing how properties such as feed composition affect the optimal filter design. In addition, we propose a novel multi-stage filtration strategy, which provides a significant mass yield improvement for the desired particles, and, surprisingly, higher purity of the product as well.

For feed containing multiple particle species, the goal may be to remove all suspended particles, but there are many applications in which the purpose of the filtration is to remove some particle species from the feed while recovering other species in the filtrate.For example, when producing vaccine by fermentation, one would want to filter the live virus out and retain the vaccine (detached protein shell of the virus, for example [23]) in the filtrate.To our best knowledge, little theoretical study has been devoted to feed containing multiple species of particles.Some experimental results are available [19,20,22], though the focus is mostly on the specific underlying application rather than mechanistic understanding of how the presence of different particle types affects the filtration process.
Thanks to recent advances in the development of fast computational tools, numerical solution of the full Navier-Stokes equations and tracking of individual particles in the feed has become a feasible approach for modeling membrane filtration [24].Several such computational fluid dynamics (CFD) studies, particularly focusing on particle deposition on the membrane, have been conducted [24][25][26][27].Such models may be very detailed, capable of tracking hundreds of millions of particles of arbitrary type and able to reproduce certain experimental data well.However, the computational demand for application-scale scenarios is extremely high; implementation of the CFD method is highly non-trivial and time consuming, and development of simpler models, which can treat different particle populations in an averaged sense, is desirable.
In earlier work [28], Sanaei and Cummings proposed a simplified model for standard blocking (adsorption of particles, much smaller than the filter pores, onto the internal pore walls), derived from first principles.The model assumes the pore is of slender shape, with pore aspect ratio defined as typical width W divided by the length D of the pore, = W/D 1, see Fig. 1.This provides the basis for an asymptotic analysis of the advectiondiffusion equation governing particle transport within the continuum framework, valid for a specific asymptotic range of particle Péclet numbers (details can be found in [28] Appendix A).The model is consistent with one proposed earlier by Iwasaki [17] based on experiments involving water filtration through sand beds, the validity of which was further confirmed in later experiments by Ison & Ives [29].
Building on that work, we recently proposed a filtration model focusing on standard blocking with quantitative tracking of particle concentration in the filtrate, which allowed for evaluation of the filtration performance of a given membrane in terms of its pore shape and particle capture characteristics, and for optimization of filtration of a homogeneous feed containing just one type of particles [30].In the present work we extend this approach to filtration with multiple species of particles in the feed.For simplicity, we consider dead-end filtration using a track-etched type of membrane.We study how the concentration ratio of the different types of particles in the feed, and the differences in membrane-particle interaction characteristics, affect the filtration process and we formulate optimization problems to determine the optimum pore shape (within a given class of shape functions) to achieve the desired objectives.To illustrate our model behavior and its application for design optimization, we explore some hypothetical scenarios of practical interest, in particular: When there are two compounds A and B in a mixture, which filter design will produce the maximum amount of purified compound B before the filter is completely fouled?Questions such as this lead naturally to constrained optimization problems: how to design a filter such that a certain large fraction of type A particles is guaranteed to be removed, while retaining the maximum yield of type B particles in the filtrate, over the filtration duration?
We propose new fast optimization methods to solve these problems, based on quantities evaluated at the beginning of the filtration, which are over 10 times faster than the method used in our earlier work [30].Motivated by some of our findings, we also propose a new multi-stage filtration protocol, which can significantly increase the mass yield per filter of the desired compound, and simultaneously improve the purity of the final product.
Many variations on the questions we address could be proposed, and the methods we present are readily adapted to a wide range of scenarios.For brevity and simplicity, however, in the present work we focus chiefly on variants of the example outlined above to illustrate our methods.For the majority of the paper we present results for the case in which the feed contains just two particle species, noting that (within the limitations of our modeling assumptions) our model is readily extended to any number of particle species (some sample results for feed containing more than two species are included in the Appendix B).
The remainder of this paper is organized as follows.We set up our two-species filtration model in §II, focusing attention on the filtration process within a representative pore of the membrane.We then outline a number of hypothetical filtration scenarios with multiple species of particles and formulate the corresponding optimization problems in §III A.
Although our optimization criteria as defined rely on simulating filtration over the entire useful lifetime of the filter, we will demonstrate the feasibility of using data from the very early stages of our simulations as a reliable predictor of later behavior, offering a much faster route to optimization, discussed in §III B. Sample optimization results will be presented in §IV.Section §V is devoted to the summary and discussion.

II. FILTRATION MODELING WITH TWO PARTICLE SPECIES
In this paper, we focus on dead-end filtration feed solution, carrying multiple different particle species, through a membrane filter.We first highlight some key modeling assumptions: we assume the particles are non-interacting (justifiable if the feed solution is sufficiently dilute); that the particles are much smaller than the pore radius; and that the pore is of slender shape, with length much larger than its width (this is the case for "track-etched" type membranes whose pores are straight and form a direct connection between upstream and downstream sides of the membrane; see, e.g., Apel [31]).We consider only one type of fouling: the so-called standard blocking mechanism, in which particles (much smaller than pores) are adsorbed on the pore wall leading to pore shrinkage; and we inherit all the additional assumptions made in deriving the standard blocking model proposed by Sanaei and Cummings [28].Under these assumptions, we set up our model for constant pressure and constant flux conditions, in §II A and §II B respectively.

A. Solute at constant driving pressure
We consider a feed solution containing two types (different physicochemical properties) of particles, type 1 and type 2, through a planar membrane filter under constant pressure.In the presentation that follows, we use uppercase fonts to denote dimensional quantities and lower case for nondimensional quantities, which will be defined in §II C when we introduce appropriate physical scalings.We assume that the membrane is composed of identical pores of circular cross-section with radius A(X, T ) (where X is distance along the pore axis), periodically repeating in a regular (e.g., square or hexagonal) lattice arrangement.Each circular pore is contained within a regular tesselating polygonal prism, which accommodates a pore of maximum radius W (0 < A ≤ W ) and height D (see Fig. 1 for example), where W D. We define the representative pore aspect ratio = W/D 1, which will be used in our particle deposition model discussed below.The incompressible feed (assumed Newtonian with viscosity µ) flows through the pore with cross-sectionally averaged axial velocity, U p (X, T ), given in terms of the pressure P (X, T ) by where K p = A 2 (X, T )/8 is the local permeability of an isolated pore (which follows from the Hagen-Poiseuille formula, see e.g., Probstein [32], consistent with our pore shape assumption).This is equivalent to a Darcy flow model with velocity U (T ) within the membrane related to U p (X, T ) via porosity Φ m = πA 2 (X, T )/(2W ) 2 , where is the membrane permeability.The flow is driven by constant pressure drop P 0 across the membrane.Conservation of mass then closes the model, giving the equation and boundary conditions governing the pressure P (X, T ) within the membrane as Extending the approach of Sanaei & Cummings [28], we propose the following fouling model equations, which assume that the two particle types are transported independently by the solvent and do not interact with each other: where C i (X, T ) is the concentration (mass per unit volume of solution) of type i particles; Λ i is a particle deposition coefficient for type i particles; and α i is an unknown (problemdependent) constant, related inversely to the density of the material that comprises type i particles.Equations ( 6) follow from a systematic asymptotic analysis (based on the small parameter defined above, see also Table II) of advection-diffusion equations for each particle species.Equation (7) assumes the rate of pore radius shrinkage (due to the particle deposition) is a linear function of the local particle concentrations at depth X, and derives from a mass-balance of the particles removed from the feed, consistent with (6).Derivations of these results for filtration of a feed with just one particle type are given in Sanaei & Cummings [28] Appendix A and [33].

B. Solute at constant flux
Here we briefly consider how the above model is modified for the same feed solution, supplied at constant flux U 0 .As fouling occurs the membrane resistance increases, hence the driving pressure must increase to maintain the same flux through the filter.Equation (2) still holds for the superficial Darcy velocity, which is now held at constant value U 0 by adjusting the driving pressure P (0, T ), giving with just one boundary condition at the membrane outlet, In this case, the incompressibility condition is satisfied automatically.Equations ( 6)-( 7) then close the model, as in the constant pressure case.

Constant pressure
We non-dimensionalize our model ( 2)-( 7) using the following scalings, with lower-case fonts indicating the dimensionless variables: where the chosen timescale is based on the deposition rate of particle type 1.The resulting non-dimensionalized equations are listed below: Eqs. ( 2)-( 5) become so that dimensionless permeability is just a 4 with the chosen scalings; and Eqs. ( 6)-( 7) take the form where λ i = 32Λ i D 2 µ/(πW 3 P 0 ) is the deposition coefficient for particle type i, ξ = C 01 /(C 01 + C 02 ) is the concentration ratio between the two types of particles, is the ratio for effective particle deposition coefficients between the two types of particles, and 0 < a 0 (x) ≤ 1 is the pore profile at initial time t = 0. Since we consider scenarios where particle type 1 is to be removed by filtration while type 2 should be retained in the filtrate, only values β ∈ (0, 1) will be considered in this paper.The model parameters are summarized in Table II for future reference.
To solve this system numerically, we first note that Eqs. ( 13)-( 15) can be solved to give Given a 0 (x), we compute u(0) via Eq.(19), which allows us to find u p (x, 0) via Eq.( 13).We then compute c 1 (x, 0), c 2 (x, 0) via Eqs.( 16) and ( 17) respectively.With c 1 , c 2 determined we then compute the pore shape a(x, t) for the next time step via Eq.( 18), then repeat the above process until the chosen termination condition (based on flux falling below some minimum threshold) for the simulation is satisfied.

Constant flux
Most scales follow from the constant pressure case of §II C 1; here we highlight only the differences for the constant flux scenario, again with lower case fonts indicating the nondimensionalized variables: The remaining scalings are as in Eqs.(11) and (12), leading to the model To solve these equations numerically, we proceed as in the constant pressure case, with the simplification that u = 1 and u p is a known function of a (22).Note that the inlet pressure p(0, t) is given by from which it follows that, as the pore radius a(x, t) decreases due to fouling, the driving pressure must increase to maintain the constant flux.We continue the simulation until the specified termination condition (based here on exhausting some fixed amount of feed, subject to a constraint on maximum inlet pressure p(0, t)) is reached.

III. OPTIMIZATION
In this section, we explore the specific scenarios introduced in §I to find the optimized initial pore shape a 0 (x) by defining a suitable objective function J(a 0 ) with corresponding constraints.For the purpose of the mathematical formulation of the optimization problem, we assume a 0 (x) ∈ C([0, 1]) (the class of real-valued functions continuous on the real interval [0, 1]); however, for practical purposes to obtain solutions within reasonable computing time we restrict the search space for the optimizer a 0 (x) to low degree polynomial functions (numerical implementation details will be given in §III B).In addition, we require 0 < a 0 (x) ≤ 1 so that the initial profile is contained within its unit prism (see Fig. 1).
In III A, we define key metrics that we use to measure the performance of the filter design and use these to set up the optimization problems.Three filtration scenarios will be studied: two under constant pressure conditions and the third under constant flux.In III B, we outline our optimization methods: first a "slow method" (described in III B 1) based directly on the objective function defined in III A below; then we propose a "fast method" (in III B 2), motivated by results obtained using the slow method.We demonstrate the feasibility of using our model with fast optimization to predict and optimize for various filtration scenarios with multiple species of particles in the feed.

A. Definitions and objective functions
Adapting the approach taken in our earlier work [30], we first define some key (dimensionless) quantities that will be used to measure the performance of the membrane.We define instantaneous flux through the membrane as u(t), and cumulative throughput j(t) as the time integral of the flux, We denote the instantaneous concentration at the outlet (x = 1) for each particle type i in the filtrate, c i (1, t) as c i,ins (t), and the accumulative concentrations of each particle type i in the filtrate, c i,acm as Let t f denote the final time of the filtration process, when the termination condition is reached.For the constant pressure case, we define this to be when the flux drops below some specified fraction ϑ of its initial value (throughout our work here ϑ = 0.1, based on common industrial practice, see e.g., van Reis & Zydney [2]); for the constant flux case, we consider t f to be the fixed time at which the specified amount of feed is exhausted, assuming that the terminal driving pressure p(0, t f ) is less than the maximum operating pressure p max for all initial pore profile functions a 0 (x) in the searching space considered.
To specify the particle removal requirement from the feed for each type of particles, we define the instantaneous particle removal ratio for type i particles, R i (t) ∈ [0, 1], as where c i,ins (t) is instantaneous concentration of particle type i at the outlet and c i (0, t) is the type i particle concentration in the feed at time t. 1 Then the initial particle removal ratio R i (0) is the fraction of type i particles removed after the feed passes through the clean filter.We also define the cumulative particle removal ratio for type i particles, Ri (t), as where c i,acm (t) is defined in (28).The final cumulative particle removal ratios at the end of the filtration are then Ri (t f ).
Preliminary investigations for our multi-species filtration model indicate that the particle removal capability of the filter improves, for the constant pressure scenarios, as the filtration proceeds and pores shrink, thus in our optimizations we impose the particle removal requirement only at the initial step, i.e., we require R 1 (0) to be greater than a specified number (R) between 0 and 1.Throughout this work we consider the desired final particle removal ratio for type 1 particles to be 0.99 and denote this fixed value by R. Other values 1 In the problems that we consider but if we wish to consider feed with time-varying particle concentrations, then c i (0, t) in ( 30) should be replaced by appropriate averaged concentrations, ci (0, t) := (

Metric Description
Range/value/definition t) instantaneous concentration at the outlet for each particle type i = c i (1, t) ∈ (0, c i (0, t)) c i,acm (t) accumulative concentrations of each particle type i in the filtrate ∈ (0, c i (0, t)) (defined in Eq.( 28)) instantaneous particle removal ratio for type i particles Table III: Key metrics defined in §III A and §IV A 2 for measuring membrane performance and their ranges, values (where fixed across all simulations) or definitions.
of R are used for "intermediate" filtration stages in our description of multi-stage filtration later, with the understanding that the final goal is to reach removal ratio R = 0.99.With these definitions, for the constant pressure case, we illustrate our methods by considering a number of membrane design optimization scenarios, outlined below.
Problem 1.In many situations there are competing demands and it may be useful to consider objective functions that assign weights to different quantities of interest.Suppose we have a feed with known concentrations of type 1 and type 2 particles, where the goal is to remove type 1 particles from the feed, while retaining type 2 particles in the filtrate and simultaneously collecting as much filtrate as possible, until the termination time t f := inf t : u(t) ≤ ϑu(0) is reached.Which filter design a 0 (x) -the initial pore profile, within our searching space -will remove a specified fraction R ∈ [0, 1] of type 1 particles and simultaneously maximize the objective function J(a 0 ) := w 1 j(t f ) + w 2 c 2acm (t f ) (where w 1 and w 2 are weights associated to the total throughput and final cumulative concentration of type 2 particles in the filtrate, respectively)?For example, in water purification [35], type 1 particles could be toxins like lead (which we insist are removed), while type 2 particles are desirable minerals.In this application, it is of interest to retain type 2 particles, but the primary concern is to produce the purified water, so a larger value might be assigned to w 1 than w 2 .This example motivates the following design optimization problem.
Optimization Problem 1 Maximize subject to Eqs. ( 15)-( 18), and Here, we seek the optimum pore shape a 0 (x) to maximize J(a 0 ), a weighted combination of j(t f ) and c 2acm (t f ), subject to the flow and fouling rules dictated by our model (Eqs.( 15)-( 18)), and the physical constraints that the pore is initially contained within the unit prism (so that adjacent pores cannot overlap), and the desired user-specified fraction R of type 1 particles is removed from the feed at the start of filtration.For example, if w 1 = 1, w 2 = 0, R = R then we are maximizing the total throughput of filtrate, with a hard constraint that at least 99% of type 1 particles are removed initially, and no concern for the proportion of type 2 particles retained in the filtrate.On the other hand, if w 1 = 0.5, w 2 = 0.5, R = R then (assuming the dimensionless quantities j(t f ) and c 2acm (t f ) are of similar magnitude) we care equally about total throughput and the proportion of type 2 particles retained in the filtrate, again with a hard constraint on removal of type 1 particles.
Problem 2. Suppose we have a large quantity of feed containing known concentrations of type 1 and type 2 particles, where the goal is to remove type 1 particles and collect the maximum quantity of type 2 particles in the filtrate (e.g., for vaccine production after fermentation, one would want to filter out the live virus -type 1 particles -and retain as much vaccine -type 2 particles -as possible in the filtrate), until the termination time t f is reached.Which filter design a 0 (x), within our searching space, will remove a specified fraction R ∈ [0, 1] of type 1 particles and simultaneously maximize the final yield of type 2 particles in the filtrate, c 2acm (t f )j(t f ) =: J(a 0 )?This question leads to the following design optimization problem.

Optimization Problem 2
Maximize subject to Eqs. ( 15)-( 18), and Here we seek the optimum a 0 (x) that maximizes objective function J(a 0 ), representing the final mass of type 2 particles in the filtrate, subject to the flow and fouling rules of our model, the physical constraints, and the desired particle removal requirement.
For the case in which constant flux through the filter is specified, we consider the following illustrative scenario: Problem 3. Given a fixed amount of feed, what is the best filter design to maximize the yield of a purified compound of interest (e.g., gold in mining [36], vaccine extraction [5] or other bio-product purification), while removing an impurity?In this scenario, we consider the optimization problem as finding the initial pore profile a 0 (x) such that the filter removes a certain proportion (R) of type 1 particles from the feed, while maximizing the amount of type 2 particles collected in the filtrate, until the feed is exhausted, i.e. the termination time t f is reached.This problem statement motivates the following optimization problem: subject to Eqs. ( 21)-( 25), and t f specified (time at which feed is exhausted).
The objective function J(a 0 ) in (33) represents choosing the optimum a 0 (x) to maximize the final mass of purified type 2 particles obtained at the end of the filtration.

B. Optimization methodology
The design optimization problems outlined above are mathematically challenging and computationally expensive in general due to non-convexity [37] (of both the objective function and the constraints), large number of design variables (in our case the number of possible design variables is infinite, as our searching space for the pore shape a 0 (x) is the infinitedimensional function class C([0, 1])) and the computational cost of evaluating the objective function (which requires that we solve the flow and transport equations until the termination time t f ) many times.For simplicity and efficiency we therefore restrict our searching space for a 0 (x) to be (low degree) polynomial functions, the coefficients of which represent our design variables (the class of searchable functions could be expanded without difficulty but with commensurate increase in computational cost).In the following two subsections we outline our two optimization routines: first the slow method, which arises naturally from the problems posed and which relies on running many simulations over the entire lifetime of the filter; then the proposed new fast method, which uses data from only the very earliest stages of filtration to predict the optimum over the filter lifetime.The two methods are compared in §IV.

Slow method
We vary the coefficients of polynomials a 0 (x) to find the values that maximize the objective functions defined in Problems 1-3, under the constraints specified in each case.
The functions a 0 (x) are referred to as shape functions in the shape optimization literature [38].In the interests of reducing computation time, for the purpose of the demonstration simulations presented here, we restrict our searching space to be the linear pore profile, i.e.
we consider initial profiles of the form a 0 (x) = d 1 x + d 0 , where d 1 , d 0 are the design variables to be optimized for each specific scenario, with searching range (d We use the MultiStart method with fmincon as local solver from the MATLAB® Global Optimization toolbox for this optimization.Since the routine finds a minimizer while we want to maximize J(a 0 ), we work with the cost function −J(a 0 ); more details of the implementation of the cost function and constraints can be found in our earlier work [30].
We specify a starting point We use the best minimizer from the list as the coefficients for our optimized linear pore profile.Note that there is no guarantee the method will find the global minimizer due to the nature of gradient descent methods applied to non-convex problems (the result found depends on the starting-points); however, local minimizers can be systematically improved and for practical purposes may be useful if they provide significant improvement over current practice (see e.g., the study by Hicks et al. on airfoil design [39]).Simulations using this method are presented in Figures 2, 3, 4, 5, described in §IV later.

Fast method
The "slow" optimization method described above is straightforward and easy to implement, but reliable results require that many (n large) individual model simulations be run through to the termination time t f .The results presented in this paper are restricted to optimizing membrane structure within the class of linear pore profiles only, but in any real application it may be desirable to optimize over wider function classes, e.g., polynomials of higher order.We find (empirically) that each unit increase of the degree of polynomial a 0 (x) requires roughly a 10-fold increase in the number of searching points to reach the best local optimum, with a corresponding increase in the run time.Run time will also increase if more than two particle species are considered, or if some of the constraints are removed or inactive (e.g., a less strict particle removal requirement) and the feasible region becomes very large.Maximum computational efficiency in practical situations is therefore critical.Motivated by the idea that imposing carefully-chosen conditions on the initial state of a system can, in many cases, guarantee certain features of later states, we propose a fast method based on simulations of the very early stages of filtration.We note that similar ideas have been used to estimate filter capacity (the total amount of feed processed during a filtration) using a method called V max , which essentially predicts the filter capacity using only the first 10-15 minutes of filtration data [40].
Extensive preliminary simulations for Problem 1, with w 2 = 0,2 indicate that the (u(t), j(t)) flux-throughput graph at optimum is initially flat and high (see Fig. 3(a) for example), with small gradient |u (0)| and large vertical intercept u(0) in comparison with graphs for sub-optimal solutions.Moreover, fouling shrinks the pore and increases resistance, thus flux decreases in time and u (0) < 0 for all model solutions.We therefore expect that, at optimum, u(0) should be as large as possible and u (0) as close to zero as possible.Similar ideas apply to the case w 2 > 0, where we wish also to maximize c 2acm (t f ), the cumulative concentration of type 2 particles in the feed at the final time: we propose instead to maximize a function based on the initial state of the system as characterized by c 2ins (0) and c 2ins (0) (respectively the initial concentration and initial concentration gradient, with respect to t, of type 2 particles at the membrane outlet x = 1).Again, preliminary simulations indicate that at optimum c 2ins (0) is close to zero and negative,3 while c 2ins (0) is large, compared to sub-optimal solutions.With these observations, we expect to maximize c 2acm (t f ) by insisting on large c 2ins (0) and small initial gradient c 2ins (0).
With these motivations, we now define modified objective functions for our fast method.
In place of (31) in Optimization Problem 1, we propose the following fast objective function, which uses data from only the initial stage of the model solution: J 1,fast (a 0 ) = w 1 u(0) + w 1 u (0) + w 2 c 2ins (0) + w 2 c 2ins (0), (34) in which the terms in w 1 act to maximize total throughput and those in w 2 maximize concentration of type 2 particles in the filtrate, where w 1 and w 2 can be tuned depending on the relative importance of the two quantities.Note that the weights assigned to u(0) (c 2ins (0)) and u (0) (c 2ins (0)) do not have to be the same; we could allow four independent weights for the four quantities in (34).However, for the simple application scenarios we considered we found just two independent weights w 1 , w 2 to be sufficient to give reliable results in an efficient manner.
To replace (32) in Optimization Problem 2, we propose the following fast objective function in which u(0)c 2ins (0) captures the initial collection of the particle 2 in the filtrate.Other forms involving u (0) and c 2ins (0) were tested, but found to confer no improvements, hence we opt for the simplest effective objective.
In the following section we demonstrate that our fast optimization method always gives results at least as good as those for the slow method, and then use it to investigate various model features and predictions.

IV. RESULTS
In this section we present our simulation results for several two-species filtration scenarios.
We focus on the effects of ξ, the concentration ratio of the two particle types in the feed, and , the ratio of the effective particle deposition coefficients for the two particle types (both these parameters are unique to multi-species filtration, having no counterparts in single-species models).For most of our simulations, we fix λ 1 = 1 (particle type 1 has fixed affinity for the membrane throughout) and the initial fraction of type 1 particles to be removed is fixed at R. In §IV A, we first present sample comparison results between the fast and slow methods for Problem 1 and Problem 2, noting that many more tests than are presented here were conducted to verify that the fast method reliably finds optima as good as or superior to those found by the slow method, under a wide range of conditions.We then use the fast method to study the effects of varying parameters β and ξ for two-species filtration under constant pressure conditions.Based on these observations, we propose a multi-stage filtration strategy that will increase the mass yield of particle we wish to recover in §IV A 2.
We present sample results for the constant flux case in §IV B, focusing on Problem 3.
A. Optimization of constant pressure filtration

Efficacy of the fast optimization method
We begin by demonstrating both slow and fast optimization methods described in §III B above. Figure 2 shows the fouling evolution of the optimized membrane pores a 0 (x) for  finds an optimized a 0 (x) as good as or better than that found by the slow method (larger or the same values for j(t f ) and c 2acm (t f ), while always satisfying the removal criterion R for particle type 1).
These results, as well as many others not discussed here, demonstrate that with the same number of searching points the fast method converges to an optimizer that in all cases is as good as or better than that obtained using the slow method, with considerably shorter running time (a typical optimization for the slow method takes 40 minutes with 10,000 searching points, while the fast method takes only 4 minutes).We also observe that varying the weights [w 1 , w 2 ] does not change the optimized profile significantly, indicating that maximizations of j(t f ) and of c 2acm (t f ) are correlated for the parameter values considered.
One possible explanation for this correlation is that, provided the type 1 particle removal constraint R 1 (0) ≥ R is met, the initial concentration of type 2 particles in the filtrate c 2 (1, 0) should be maximized by maximizing the initial flux u(0), since the higher the flux, the more type 2 particles will escape capture by the filter.
Figure 4 presents direct comparisons of the slow and fast methods for Problem 2. Results for the slow method, with objective function (32), are indicated by dashed curves; and those for the corresponding fast method, with objective function (35), by dotted curves.
The left panel, Figs. 4 (a-c), shows results for various feed particle-composition ratios ξ (other parameters as in Fig. 3); while the right panel, Figs. 4 (d-f), compares results for various effective particle-membrane interaction ratios β, with ξ = 0.5.The flux through the membrane and the cumulative particle concentrations of each particle type in the filtrate are plotted as functions of filtrate throughput over the duration of the filtration, 0 ≤ t ≤ t f .In all cases, for the same number of searching points, the fast method converges to the same optimal pore profile as the slow method across all ξ and β values considered (though the optima obtained are different for each parameter set).Similar to Problem 1, the computational speedup is considerable using the fast method.
In addition to demonstrating the efficacy of the fast optimization method, the results also illustrate some general features of the model.When a feed contains a larger fraction (higher ξ-value) of particles to be removed (type 1 particles here), our model predicts shorter filter lifetime (due to faster fouling) when compared to a feed with lower ξ-value, leading to less total throughput and lower final accumulative particle concentration of type 2 particles in the filtrate, see e.g., Figs.4(a) and (c).This is not desirable if we want to maximize total collection of type 2 particles; we will present one possible way to circumvent this issue in §IV A 2, where a multi-stage filtration is proposed.
Figure 4(b) is the (c 1acm , j) plot.Note that in all cases, the constraint for removal of particle type 1 is tight at the optimum, with the exact specified proportion R of particles removed from the feed at time t = 0. Figures 4 (d-f) show that larger β values (meaning that the two particle types are more physicochemically similar; recall β ∈ (0, 1) throughout our study, and if β = 1 both particle types interact identically with the membrane) lead to faster fouling of the filter, with lower total throughput and lower total yield of type 2 particles in the filtrate.This confirms our expectation that the more similar the particle types are, the more challenging it is to separate them by filtration.To achieve effective separation, a sufficient physicochemical difference γ = 1 − β is required.
Encouraged by the excellent results and significant speedup obtained when using the fast method with the same number of searching points as the slow method, we next investigate its performance with fewer searching points.Figure 5 shows the comparison between the slow method (dashed curves) with 10,000 searching points (found, empirically, to be the minimum number required for reliable results) and the fast method (dotted curves) with decreasing number of search points (10,000, 1,000, 100).Model parameters are fixed at ξ = 0.5 and β = 0.1; other parameters are as in Fig. 3.These results (as well as many other tests, not shown here) indicate that the fast method produces reliable results with just 1,000 searching points (blue dotted curves; this optimum even appears superior to the slow method with 10,000 search points, providing slightly higher total throughput).Even with as few as 100 search-points the fast method produces reasonable, though suboptimal, results (black dotted curve).In all cases the particle removal constraint on c 1 is again tight at the optima found.Since run time for the optimization routine scales in direct proportion to the number of searching points, a 10-fold reduction in the number of search points needed represents a significant additional computational saving: the fast method with 1,000 search points is approximately 100 times faster than the slow method with 10,000 points.

Multi-stage filtrations
In Fig. 4 (a) we observed that, with a higher concentration ratio of type 1 particles in the feed, the optimized filter for a single-stage filtration tends to be fouled faster, which leads to lower total throughput per filter.This makes sense as the filter needs to remove a higher mass of impurity (type 1 particles) to achieve the initial particle removal threshold R 1 (0) ≥ R when ξ is larger.Our simulations also reveal that the fouling is largely confined to a narrow region adjacent to the upstream surface of the filter at optimum, with the majority of the downstream portion of the filter remaining unused (see Fig. 2).In this section we propose a multi-stage filtration scenario that has the potential to alleviate some of these inefficiencies.
Heuristically, we would like to process more feed per filter by increasing the membrane porosity and simultaneously make more efficient use of the membrane material by fouling a higher proportion of the pore (void) volume.However, increasing the porosity in general decreases the particle removal efficiency, so the filtrate will require further purification to satisfy the particle removal requirement.This provides the motivation for the proposed multi-stage filtration strategy: we will lower the initial particle removal requirement to increase the amount of feed processed per filter and try to satisfy the particle removal requirement by filtering the collected filtrate again, possibly more than once (multi-stage).
The multi-stage filtration will be cost-effective if the increase in feed processed can offset the increase in the number of additional filters required to meet the particle removal requirement.
From the optimization point of view, we increase the feasible searching space by relaxing the initial particle removal constraint so that a better optimizer might be found.
In the following discussion we focus on the optimization Problem 2, where the goal is to maximize the yield of type 2 particles per filter used, while achieving effective separation, 4which for definiteness we here define as removing the desired fraction R of particle type 1 from the feed ( R1 (t f ) ≥ R) while simultaneously recovering a minimum desired yield fraction Υ of type 2 particles in the filtrate ( R2 (t f ) ≤ 1 − Υ).For all of our simulations, R = 0.99 and Υ = 0.5.We define the purity for type i particles in the filtrate, k i ∈ [0, 1], as where c i,acm (t f ) is the accumulative concentration of the type i particle in the filtrate at the end of the filtration.With our hypothesized scenario of feed containing desired (type 2) and undesired (type 1) particles in mind, we note a simple relationship between the purity of type 2 particles and the final cumulative removal ratios: We will return to these definitions in our discussion of the multi-stage filtration results below.
The basic idea behind our multi-stage filtration is to first optimize the filter with a less strict initial type 1 particle removal requirement (i.e., we require R 1 (0) ≥ R < R) and filter the feed solution two or more times to achieve a larger total yield per filter of purified type 2 particles than in a single-stage filtration, with the effective separation condition satisfied at the end of the multi-stage filtration.We determine the stage of filtration by how many times the solution has passed through clean filters: for example, the clean stage 1 filter will take feed directly and be used to exhaustion; the filtrate collected from the stage 1 filter will then be sent through a new (clean) stage 2 filter (which may be used more than once within stage 2).
We propose the following two-stage or multi-stage filtration strategy for Problem 2 (also summarized as a flow chart in Figure 6): 1. Optimize the filter with a less strict initial particle removal requirement than desired (R 1 (0) ≥ R < R); denote the optimized filter as F R (e.g., for R = 0.5, we denote the optimized filter as F 0.5 ). 2. (Stage 1) Run the filtration simulation using F R until the filter is completely fouled; collect the filtrate.In scenarios to be considered later we allow stage 1 to use several filters simultaneously.In order to keep track of the number of filters used and how many times each is reused, we identify each (stage) filter used by F R,m , e.g., the second (stage) F 0.5 filter will be denoted , and we track how many times each stage filter has been used by n(m), e.g., n( 2) denotes the number of times F 0.5,2 is used.In cases where stage 1 involves more than one clean filter (used simultaneously) we also use notation l m to denote the total number of mth stage filters used, e.g., l 2 is the number of stage 2 filters used; and we denote the total number of filters used at the end of multi-stage filtration by M = m l m .After the multistage filtration is concluded we calculate the total mass of type 2 particles collected per filter used, c 2acm (t f )j(t f )/M , and compare it with the collected mass from the filter F R optimized for single-stage filtration.If the mass collected per filter used is larger for the multi-stage filtration, then the process is deemed more cost-effective.
In table IV, we list results comparing a single-stage filtration using a filter F R (optimized for particle type 1 initial removal threshold set at the desired value R), with two separate two-stage filtrations using filters F 0.7 and F 0.5 (optimized for lower thresholds R = 0.7 and R = 0.5, respectively).In all cases fast optimization was carried out using objective function J 2,fast (a 0 ) = u(0)c 2 (0) with ξ = 0.9, β = 0.1, α 1 = α 2 , and λ 1 = 1.The quantities listed in table IV are: R, initial type 1 particle removal threshold; M , total number of filters used in each case; n(2), the number of times the 2nd (stage 2) filter is used; R1 (t f ) and R2 (t f ), the final cumulative particle removal ratios for particle types 1 and 2 respectively; k 2 , the purity of type 2 particles in the final collected filtrate; j(t f ), the total throughput; and the total mass yield of purified type 2 particles per filter (all relevant quantities are defined in Table III).Table IV: Comparisons of single-stage filtration (R = R) with 2-stage filtrations (R = 0.7 and R = 0.5).
We record: M , the total number of filters used for each filtration process; n(2), the number of times the 2nd filter is used for each multi-stage filtration process; R1 (t f ) and R2 (t f ), the final cumulative particle removal ratios for particle types 1 and 2 respectively; k 2 , the purity of type 2 particles in the final filtrate; j(t f ), total throughput; and type 2 particle mass yield per filter.
These preliminary results show that, when our multi-stage filtration protocol is applied, we can achieve the same final particle removal requirement R1 (t f ) ≥ R as the single stage filtration but with much higher yield per filter of particle type 2; the third row of Table IV shows that, with R = 0.5 the yield of purified type 2 particles per filter is almost doubled when compared to the single-stage filtration optimized for R = R. From table IV we also observe that the multi-stage filtrations improve the purity of the filtrate as indicated by the k 2 values.We note that all three filtrations achieve effective separation according to our (somewhat arbitrary) definition, which corresponds to purity k 2 ≥ 0.847 for the cases considered in table IV, i.e. ξ = 0.9.If higher purity is desired to consider a separation effective, the removal ratios R1 (t f ) and R2 (t f ) can be adjusted accordingly based on Eq. (37).
Figure 7 illustrates the results for the optimized filters summarized in table IV and described above via the fouling evolution of the filter pores.Figures 7 (a-c filters, F 0.7,2 and F 0.5,2 respectively.We can see that when the initial removal threshold R is decreased, the fouling of the pore becomes more uniform along its depth, and the porosity of the corresponding optimized filter F R increases.For the case R = 0.5, the optimized pore profile is almost as wide as possible; the gray colored region corresponding to the membrane material is too thin to be visible.The high mass yield per filter and small quantity of membrane material required to produce F 0.5 indicates that if the membrane material has good selectivity (higher λ 1 value in our model), it might be advantageous to focus on maximizing filter porosity as a design approach to increasing the mass yield per filter, while simultaneously reducing the membrane material cost per filter and achieving effective separation using multi-stage filtration.
From Figs. 7 (d) and (e) it is clear that in both multi-stage filtration protocols, the secondary filters F 0.7,2 and F 0.5,2 are only lightly-used at termination, and could be used to process more filtrate.Specifically, we could use two or more first-stage filters F R,1 in order to create sufficient once-filtered fluid to pass through the second stage filter F R,2 and foul it significantly.We anticipate that increasing the volume of filtrate collected from stage 1 of the filtration using multiple F R,1 should lead to higher mass yield per filter by more fully utilizing the filtration capacity of the stage 2 filter F R,2 .Before investigating this idea in detail we first test it using two, three and four stage-1 filters (l 1 = 2, 3, 4), which lead to 2-stage, 2-stage and 3-stage filtrations, respectively.The results are presented in Fig. 8, using the filter F 0.5 optimized for R = 0.5 as in Fig. 7.
In the first test example, in stage 1 we collect filtrate by exhausting two F 0.5,1 filters (l 1 = 2; the fouling plot for each of these F 0.5,1 is identical to Fig. 7 (c) so is omitted in Fig. 8) and then in stage 2, we send the combined filtrate repeatedly through an initially clean F 0.5,2 .Figure 8(a) shows the subsequent fouling of F 0.5,2 , with alternating blue and red color indicating particle deposition and filter reuse as before.After passing the filtrate through F 0.5,2 three times, the final particle 1 removal requirement is met (so l 2 = 1 and In the second example, in stage 1 we collect filtrate by exhausting three F 0.5,1 filters (l 1 = 3; again see Fig. 7(c) for this stage).In stage 2 we pass the filtrate from stage 1 through an initially clean F 0.5,2 .Fig. 8(b) shows the fouling of this second-stage filter F 0.5,2 .
It is used twice, but during the second use becomes completely fouled before all the filtrate can be filtered.Leaving aside for the moment the question of whether a second stage-2 filter should be introduced to deal with the leftover twice-filtered fluid, we check the (thrice filtered) filtrate from this second stage-2 filtration and find that it meets the final particle 1 removal requirement.In this example, l 2 = 1 and M = 1 1 + l 2 = 4.
In the third example, at stage 1 we collect filtrate by exhausting four F 0.5,1 filters (l 1 = 4).
This combined filtrate is then passed through a clean second-stage filter, F 0.5,2 , whose fouling is shown in Fig. 8(c).This F 0.5,2 filter is completely fouled after one use (l 2 = 1).Again, we defer the question of whether a second stage-2 filter would be cost-effective to deal with the remaining once-filtered feed, and check the particle-1 removal requirement of the twicefiltered feed.It is not yet satisfied, so we need a third stage of filtration with a new filter F 0.5,3 .Fig. 8(d) shows the fouling of this F 0.5,3 filter, which is used three times before the final particle 1 removal requirement is satisfied (l 3 = 1).Here M = l 1 + l 2 + l 3 = 6.
We find that in the first example, when we collect filtrate from two F 0.5,1 filters (Fig. 8(a)), the mass yield of type 2 particles per filter is 0.012, which is indeed greater than the value 0.0091 obtained with the original two-stage filtration of Fig. 7.However, with three stage-1 F 0.5,1 filters, the second example of Fig. 8(b), the mass yield of type 2 particles per filter decreases to 0.0067 (see table V), which may be explained by the fact that the second stage filter F 0.5,2 is completely fouled on its second use before all the filtrate obtained in stage 1 can be processed (the yield loss is due to the discarded filtrate).Similar loss of filtrate is observed in the third example, the 3-stage filtration of Figs.8(c) and (d), in which filtrate collected from four F 0.5,1 stage 1 filters was sent through a stage-2 filter F 0.5,2 , which is exhausted before all of the stage-1 filtrate can be filtered a second time.Despite this loss, the mass yield per filter is 0.010, nearly as good as the first example of Fig. 8(a).Additional simulations of the second and third test scenarios, in which new filters were introduced to process the discarded filtrate, gave less favorable results than those presented here.These three multi-stage filtration experiments suggest that a single stage-2 F 0.5,2 filter can process filtrate collected from three to four stage-1 F 0.5,1 filters, but no more.
The observations of Fig. 8, though preliminary, indicate there may be an optimal ratio between the number of filters to use at different stages, which would utilize each filter's filtration capacity as fully as possible, and minimize the loss of filtrate at each stage, ultimately maximizing the mass yield per filter.We used our model to conduct such an investigation, the details of which are provided in Appendix §A.We find (by trial and error) that the mass yield per filter can be as high as 0.013 using a four-stage filtration, with the following numbers of filters per stage: l 1 = 18, l 2 = 6, l 2 = 3, l 4 = 1.This four-stage filtration is illustrated schematically in figure 9.

B. Optimization of constant flux filtration
In this section we briefly highlight results for the constant flux case, focusing on optimization Problem 3, with objective function (33).We study how the particle composition ratio in the feed (ξ) and the quantity of feed processed affect results.In all simulations we impose the additional constraints that the initial driving pressure p(0, 0) ≤ 100 and the driving pressure at the end of the filtration should be no greater than 10 times the initial driving pressure, i.e., p(0, t f ) ≤ 10p(0, 0) (typical driving pressure at termination is about 1.5-2 times the initial pressure [44,45]).Simulations of our model at constant flux show that particle retention typically deteriorates over time, therefore, in addition to the initial particle removal requirement (R 1 (0) ≥ R) we also impose that the accumulative particle removal should be greater than a fixed number at the end of the filtration ( R1 (t f ) ≥ R, where R ≤ R; this requirement also means that no fast optimization method is practicable for this case).The quantity of feed processed is fixed by specifying the number of time iterations N (in all simulations presented here N ≤ 1000 and the time step is fixed).
Figure 10 shows the evolution of the optimized pore profiles obtained with ξ = 0.9, 0. profile is of Λ shape (instead of the V shape we observed consistently in the constant pressure case) and the particle deposition is more evenly distributed over the length of the pore.
Also, in contrast to the constant pressure case, we see that the particle type 1 retention capability of the filter decreases in time for all three ξ values, with the most significant deterioration observed for the feed containing the highest fraction of impurity (the largest ξ-value, ξ = 0.9, see Fig. 10(b)).
Another question of interest for this constant flux case is: how does the amount of feed processed affect the optimization result?We illustrate this by considering three different values of N , the total number of timesteps in our simulations.Figure 11 shows the evolution of pore profiles optimized for the constant flux objective function (33), with ξ = 0.9, β = that, as the quantity of feed decreases, the optimized pore profile changes from a V to a Λ shape.Comparing the optimized pore profiles with Figs.11(a-c) we see that the Λ shape is more prone to driving pressure increase and particle retention deterioration, as well as the more even distribution of fouling noted earlier; observations that we now explain.
We deal first with the observation that for pores of Λ shape, particles deposit more evenly along the pore depth compared to pores of V shape.Particle concentration is always highest at the pore entrance, which favors a high deposition rate; however, flux u p is also highest here for pores of Λ shape, which is unfavorable for particle deposition (both observations follow from Eq. ( 23)).On the other hand, at the pore exit, particle concentration is lowest (unfavorable for deposition); but flux u p is also lowest (favorable for particle deposition).
Hence, for pores of Λ shape, there is always a competition between particle concentration and flux, which leads to the observed even fouling distribution along the pore length.
We next argue heuristically that this more uniform particle deposition is responsible for the observed particle concentration increase as follows: From Eqs. ( 22), ( 23) and ( 24) we obtain showing that the change in particle concentration at outlet for type i particles depends on the change in the value of 1 0 a(x, t)dx.For Λ-shaped pores particle deposition is more even, thus a(x, t) changes over the entire depth of the pore, with the consequence that 1 0 a(x, t)dx changes more significantly than for V-shaped pores, where a(x, t) changes significantly near the pore entrance, but on a region of small measure.The net effect for the Λ-shaped pore is the observed particle concentration increase in time in the filtrate.The same argument may also explain the significant pressure change for pores of Λ shape compared with pores of V shape, as the pressure change depends on the change of 1 0 a −4 (x, t)dx, see Eq. (26).Collectively, these arguments suggest the following explanation for why the Λ shape is selected for lower quantities of feed: For less feed, the filtration duration will be shorter; the significant particle concentration increase at the beginning of filtration observed with the Λ-shaped pore (see Fig. 11 (c)) is favorable for increasing the mass yield of type 2 particles, while the short filtration duration keeps the concentration increase for type 1 particles within the prescribed removal limit.
The results in Figure 11 raise the question of whether the optimized pore profile will take the Λ shape more generally for sufficiently small feed quantity at constant flux regardless of particle composition ratio ξ. Figure 12 shows a sequence of simulations with a small quantity of feed characterized by N = 100, for different feed particle composition ratios ξ = 0.9, 0.5, 0.1, with β = 0.1, λ 1 = 10, R = 0.98.The figure shows: (a) (p(0, t), j)-plot; (b) (c 1acm , j)-plot; (c) (c 2acm , j)-plot; while (d,e) show the pore profiles at the termination of the filtration, with blue color indicating deposited particles, for (d) ξ = 0.5; (e) ξ = 0.1 (the corresponding result for ξ = 0.9 was shown in Fig. 11(e)).Collectively, Fig. 11(e) and Figs.12(d) and (e) suggest that, for sufficiently small feed quantity, the optimized pore profile takes a Λ-shape regardless of feed particle composition.From Figs. 12(b) and (c) we see the particle concentration changes are not significant for these short duration cases (N = 100), and the particle removal requirement for type 1 particles is satisfied for all three feeds with different particle-composition ratios; however, the pressure increase is still visible, see Fig. 12(a).

V. CONCLUSIONS AND FUTURE STUDY
In this work we proposed a simplified mathematical model for filtration of feed containing multiple species of particles.Our focus in the main body of the paper was on a feed that contains just two particle species; a brief discussion of how the model extends to an arbitrary number of species is given in Appendix B. For the two-species case, two important model parameters were identified and investigated to elucidate their effect on separation and optimal filter design: ξ, the concentration ratio of the two particle types in the feed, and β = Λ 2 α 2 /(Λ 1 α 1 ), the ratio of the effective particle deposition coefficients for the two particle types.A number of optimization problems for maximizing the mass yield of one particle species in the feed, while effectively removing the other, were considered, under both constant pressure and constant flux driving conditions.For filtration driven by a constant pressure drop, we found that the optimized pore profile is always of V-shape, which is in agreement with our earlier findings [30] for single-particle-species filtration (where the goal is to maximize total throughput of filtrate over the filter lifetime while removing a sufficient fraction of impurity).For filtration driven by a constant flux, the optimized pore profile may take either a V-shape or a Λ-shape depending on the particle composition ratio and the amount of feed considered for the optimization scenarios.
To increase the appeal and utility of our model for filter design applications, we proposed new objective functions (the fast optimization method) based on evaluating key quantities at the initial stage of the filtration.Due to the simpler forms of the proposed objectives, the fast method can be carried out with a relatively small number of initial search-points in design parameter space (compared with the slow method, which requires that a large number of simulations be run through to filter failure time).The proposed fast method is approximately 100 times faster than the naive slow method.The ideas that motivated our fast method could potentially be usefully applied to other optimization problems that require evaluation of quantities at the end of the time evolution, provided those quantities exhibit some monotonicity over time.
Observing that (based on our model predictions), effective separation in a single-stage filtration is usually achieved at the expense of short filter lifetime and inefficient filter use (most of the filter remaining only very lightly fouled), we also proposed an alternative approach for maximizing the mass yield per filter while achieving effective separation, using multi-stage filtration.With this approach we found that the mass yield per filter could be as much as two-and-a-half times that produced by the optimal single stage filtration, and surprisingly the purity of the final product is higher as well.In addition to the higher mass yield, the filter optimized for multi-stage filtration also requires less material to manufacture, due to its higher porosity.Multi-stage filtration has been utilized in industry [46] and reported experimentally [20,47]; however, to our best knowledge, little attention has been paid to optimizing this process from the theoretical side.We hope that our work will inspire further systematic studies into this promising approach.l 1 l 2 l 3 l 4 c 1acm (t f ) c 2acm (t f ) j(t f ) mass yield/filter 1 1 (4)  from the National Science Foundation under grant NSF-DMS-1615719.
Figure 13 shows, for each multi-stage filtration considered, the total mass yield per filter (Fig. 13(a)); the final purity k 2 of the filtrate (Fig. 13(b)); and the final cumulative particle removal ratios R1 (t f ), R2 (t f ) for the two particle types.We find that the maximal mass yield per filter is 0.013; the corresponding (four-stage) filtration is apparent as the global maximum of the mass yield per filter in Fig. 13(a).This maximum yield is also indicated in red font in table V, and is almost two and half times the yield obtained with a single-stage filtration optimized to maximize yield while immediately satisfying the purity constraint.
This four-stage filtration is illustrated schematically in Fig. 9.We note that the higher mass yield per filter appears to be achieved at the expense of lowest purity k 2 = 0.873 among other multi-stage filtrations considered, with final cumulative particle 1 removal ratio R1 (t f ) just above 0.99 (see Fig. 13(b,c)), which makes sense as optimizers are generally found at the boundary of the feasible search space where one or more constraints are tight.However, the 2-stage local maximum simultaneously achieves high mass yield per filter together with the highest purity, k 2 = 0.997, which indicates that this two-stage filtration may be useful to achieve high mass yield without sacrificing the purity of the final product.

Appendix B: Multiple species
In this section, we present some sample results for feed containing more than 2 species of particles.We non-dimensionalize our model ( 2)-( 7) using the same scalings as §II C 1 for most quantities, with the following variations: Eqs. ( 13)-( 15) remain unchanged, and Eqs. ( 6)-( 7) take the form where λ i = 32Λ i D 2 µ/(πW 3 P 0 ) is the deposition coefficient for particle type i, ξ i = C 0i / i C 0i is the concentration ratio of type i particles, β i = Λ i α i /(Λ 1 α 1 ) is the effective particle deposition coefficient for particle type i (relative to particle type 1), and a 0 (x) is the pore profile at initial time t = 0. To illustrate the model we consider an optimization problem similar to Problem 2. For definiteness, we consider a feed containing three species Table VI: Three species feed filtration.We record ξ i , particle ratios for each type of particles, Ri (t f ) the final cumulative particle removal ratio for particle type i, k 2 purity of type 2 particles in the final filtrate, j(t f ) total throughput.

Figure 1 :
Figure 1: Sketch of a cylindrical pore of radius A(X, T ) and length D inside a square prism, representing a basic building-block of the filter membrane (our model is relevant for any other regular tesselating pore-containing prism, e.g. a hexagonal or triangular prism).Blue arrows indicate the flow direction; colored dots indicate the different particle types present in the feed, and W represents the maximum possible pore radius.
cumulative particle removal ratio for type i particles ∈ [0, 1] R desired final cumulative particle removal ratio for type 1 particles 0.99 Υ desired fraction of type 2 particles in filtrate (effective separation) 0.5 ϑ flux fraction at termination (constant pressure filtration) 0.1 k i purity for type i particles in the filtrate at the end of filtration ∈ [0, 1] γ effective physicochemical difference between the two species ∈ [0, 1] guess for running the local solver fmincon), cost function (based on our objective functions and constraints), design variable searching range, and number of searching points n (the number of points in (d 1 , d 0 )space that will be explored) for the MultiStart method.With the user-specifed starting point (d 0 1 , d 0 0 ), an additional (n − 1) starting points (d i 1 , d i 0 ) ∈ [−1, 1] × [0, 1], i = 1, 2, ...n − 1 are generated by the MultiStart algorithm.The resulting n points are then used to run the local solver fmincon (based on a gradient descent method) to find a list of local minimizers.

Problem 1
with w 1 = 1, w 2 = 0, ξ = 0.5, β = 0.1, α 1 = α 2 , λ 1 = 1, and R 1 (0) ≥ R, using the slow method (left panel) and fast method (right panel).The top row shows the clean, unfouled optimized pore profiles at t = 0 (Figs.2 (a) and (d)); the center row shows the fouling of these pores at t = t f /2 halfway through the filtration (Figs.2(b) and (e)); and the bottom row shows the fouled pores at termination time t = t f (Figs.2 (c) and (f)).The gray region is the filter material, and the dark blue color indicates the fouling by deposited particles.The white area denotes the open pore (void), and the red center line is the axis of symmetry of the pore (which has circular cross-section).This figure illustrates that the optima a 0 (x) found by fast and slow methods are indistinguishable.

F
R used, to compute the mass yield per filter.
) show the filters from the first filtration stage, optimized for particle removal thresholds R = R (a), R = 0.7 (b), and R = 0.5 (c), at time t = t f (when the flux is reduced to the fraction ϑ = 0.1 of its initial value).The blue and red colors indicate deposited particles; a change of color indicates reuse of the filter.Figures 7 (d) and (e) show the fouling of the second stage Figure 10 shows the evolution of the optimized pore profiles obtained with ξ = 0.9, 0.5, 0.1, β = 0.1, λ 1 = 10, R = 0.98, N = 1000.The figure shows: (a) driving pressure vs throughput, (p(0, t), j)-plot; (b) accumulative type 1 particle concentration in the filtrate vs throughput, (c 1acm , j)-plot; (c) accumulative type 2 particle concentration in the filtrate vs throughput, (c 2acm , j)-plot; and (d-f) show the pore profiles at the termination of the filtration for the three ξ-values, with blue color indicating deposited particles: (d) ξ = 0.9; (e) ξ = 0.5; (f) ξ = 0.1.We observe that, as ξ varies, the optimized pore profile changes significantly.For feed containing less impurity (the smallest value, ξ = 0.1, Figure 10(f)) the optimized pore

Figure 13 :
Figure 13: (a) Mass yield per filter, (b) purity of type 2 particles (k 2 ), (c,d) final cumulative particle removal ratios for particle types 1 and 2 ( R1 (t f ) and R2 (t f )) are plotted against the total number of filters used in each multi-stage filtration.The local maximum mass yields per filter for 2-stage, 3-stage and 4-stage filtrations are labelled with a list of values (l (k) m ), representing the number of filters l m used for stage m, and the number of times k each filter is used, listed in order of increasing m.

Table II :
Dimensionless parameters and descriptions (from TableI).
3. Re-filter the collected filtrate through another clean F R ; collect the new filtrate.4.Test the filtrate from step 3. Does it meet required type 1 particle final removal requirement R1 (t f ) ≥ R?If yes, we are done; if no, repeat step 3 using the same filter until the requirement is met, or until this F R is completely fouled.5. (multi-stage) If F R is completely fouled before R1 (t f ) ≥ R, use another clean F R and repeat step 3 until R1 (t f ) ≥ R. 6.Once the threshold R1 (t f ) ≥ R is met, record the total mass of compound 2 in the filtrate and the number of R1 (t f ) R2 (t f ) k 2 j(t f ) mass yield/filter

Table V :
Comparisons of multi-stage filtrations (up to four stages are considered) with differing ratios of the number (l m ) of filters F 0.5,m used at stage m.The first four columns list values l m , with superscripts (k) indicating that each filter is used k times.The remaining columns show final cumulative particle concentration for type 1 and type 2 particles, c 1acm (t f ) and c 2acm (t f ), total throughput j(t f ) and compound 2 mass yield per filter.The global maximum mass yield per filter is highlighted in red font.