146 th European Study Group with Industry/ Co-creation event with society Optimizing the performance of a conical ceramic membrane

Smart Separations Limited (SSL) is a UK-based start-up who have developed a ceramic membrane with micron-sized conical pores distinct to the cylindrical pores typically used for ﬁltration. This new technology has the potential to be highly beneﬁcial to many applications. However to realize its potential, a comprehensive analysis of the performance and eﬃciency of the membrane is vital. We use mathematical modelling to explore and quantify the behaviour and performance of the membrane and its link to the underlying pore structure. We derive a reduced model based on the slenderness of the membrane


Introduction
Membrane separation is a vast industry with a wide range of applications, including water treatment [16,17], biopharmaceuticals [3,15] and food processing [6,8]. For example, filters are crucial to remove waste and excess water from the blood in kidney dialysis [11], and yeast and bacteria in beer production [7]. Despite the diverse industrial applications, an overarching goal of membrane design is to maximize the product yield. The process of filtering particles with diameters between 0.1-100 µm is called microfiltration. This is a poorly explored field with many important applications from the filtering of blood to the purification of water and air. Mathematical modelling can offer key insight into the filtration process and operating conditions, and thus provide a cost-effective way to optimize filter design.
A common form of membrane separation is depth filtration, in which small particles (contaminants) are trapped within, and not just on the surface of, the porous filter material. Such filters often capture the majority of the contaminants in the initial portion of the filter while leaving the latter portions relatively unused, leading to premature clogging and reduced filtration efficiency [5].
SSL fabricate ceramic membranes with micron-sized conical pores that offer great potential in the field of separation of different media/particles, with applications for the filtration of beers and wines, blood, and removal of dust and pollen from the air. The membranes are 3-5 mm thick, and can have pore radii in the range 8-100 µm. A fundamental challenge is how to design the pore structure to maximize performance in a given filtration scenario. For instance, pores that constrict with depth can offer superior removal performance. However, this must take into account design constraints. For instance, a critical maximum porosity and a minimal overall thickness are enforced to ensure the filter does not collapse during operation, while pores with larger angle are less likely to survive during the manufacturing process. The balance of energy minimisation and design constraints therefore suggests there may be an optimum pore angle and number of pores.
Mathematical and computational methods are very useful for investigating the filtration process at a fraction of the full experimental cost. Full computational fluid dynamics software provides comprehensive insight, but is time consuming. Further, although filtration occurs on the scale of the particle or pore size, it is generally the overall macroscale behaviour, such as the total mass of particles removed, which is the main concern in filter design. As such, reduced continuum models provide a balance between computational expense and insight (see, for example, [4]). In this report we shall derive a series of reduced models of this nature that provide predictions for the filter behaviour while allowing for efficient parameter sweeps to determine optimal operating regimes.
In Section 2 we outline the challenges presented by SSL at the 146th European Study Group with Industry(ESGI146). In Section 3 we present the underpinning mathematical model for the filtration process. We then use a combination of analytical and numerical techniques to solve these in a series of distinct yet complementary regimes. Specifically, in Section 4 we exploit the fact that the pores are slender to derive a reduced form of the model. In Section 5 we relax this constraint and solve the full 2D system using finite-element methods. In Section 6 we propose a new probabilistic model based on expectation values for predicting the flux decline with time due to blocking. In Sections 4-6 we solve the models to predict the behaviour of the filtration device and determine the operating regimes that lead to efficient filtration. Finally, in Section 7 we conclude our findings in a concise manner and outline the proposed next steps.

Description of the challenges
SSL are interested in understanding the behaviour and performance of their filters. In particular, they wish to explore the following:

Flow rate prediction and maximisation
What is the relationship between the air or liquid flow rate through the filter and the pore spacing, top pore radius, and conical pore angle? What is the critical concentration of pores of a specific diameter per unit area and the contribution of the remaining porosity of the skeleton?
• Given a filter structure (characterised as in Challenge 1), how fast will the membrane become clogged with different types of particles? What role does the fluid being filtered play (e.g., viscosity, temperature)?
• What role does volume of required fluid play? e.g., higher volumes of beverages are required than biological fluids.
• How can we modify the filter protocol (e.g., changing membrane orientation relative to the flow field) to reduce clogging?
• How do the filters behave when operating in a diffusive rather than convective mode? For example, how do glucose and H 3 O + molecules diffuse between two chambers separated by our membrane filters? This mode may also be used to selectively filter elements (e.g., bacteria, red cells, white cells, carcinogenic cells).

Mathematical model
The pores are of conical shape, with radius at the inlet r in and radius at the outlet r out . We define the geometry of the pores in cylindrical coordinates (r, θ, z) in three dimensions, with r = 0 corresponding to the centreline of the pore, z = 0 at the inlet and z = L at the outlet of the pore ( figure 1). An incompressible and Newtonian fluid of viscosity µ enters the pore with velocity u. In what follows, we assume that the pores and the fluid flow are axisymmetric and we formulate a two-dimensional mathematical model in the computational domain Ω (see figure 1). The fluid flow is described by the incompressible Stokes equations Here, u(r, z, t) = ue r + we z denotes the fluid velocity, where e r and e z are unit vectors in the r and z directions respectively, and p(r, z, t) denotes the pressure. The operator ∇ is defined in the axisymmetric coordinates. In addition to the Stokes equations above, we impose the following boundary conditions where equation (2a) expresses the symmetry boundary condition on the centreline and equation (2b) express no slip and no flux at the pore walls. The flow rate into the pore is defined by ESGI146/co-creation event with society (ESGI146) The fluid is contaminated with particles whose transport is described by the advection-diffusion equation where c is the contaminant concentration and D is the diffusion coefficient. We also supplement the transport equation with the following boundary conditions at the centreline of the pore and the wall where K 1 is the adsorption rate of contaminants onto the wall via diffusive transport and n is the normal to the wall. The boundary condition (5b) describes the reaction process according to which particles adhere to the wall at a rate proportional to their local concentration. The deposition of particles at the wall changes the shape of the pore wall at r = h(z, t), according to the equation where K 2 is the adsorption rate that defines how the adsorbed contaminant changes the pore geometry. This equation assumes that the pore radius shrinks in time in response to the deposition, with linear rate K 2 . The problem is written in dimensionless form by performing the following transformation of variables with W the axial fluid velocity at the inlet, R = (r in +r out )/2 the average pore radius, and C the particle concentration at the inlet. The pressure scale is chosen as P = µWL/R 2 in order to balance the axial component of the momentum equation (1a). For the time scale T , there is more than one choice but here we choose the diffusive time scale T = L 2 /D. The tildes are suppressed henceforth and the following length and velocity ratios are introduced: where the latter condition is such that to balance the continuity equation (1b). We will later solve the governing equations for two different cases: (a) in the long-andthin approximation, in which case the ratio ϵ 1, and (b) in the general scenario with ϵ = O(1).
The dimensionless model is given below in equations (9)-(11): the fluid velocities and pressure are described by with the boundary conditions (2). Here, subscripts denote differentiation. The particle concentration satisfies and the pore wall deformation is with dimensionless parameters corresponding respectively to: the Péclet number, which measures the relative rate of advective to diffusive transport, an adsorption parameter corresponding to the strength of removal of contaminant, and a shape change parameter corresponding to the effect of the particle deposition on the shape of the pore. In this report we are interested in initial pore profiles that are conical in nature. We thus define the initial profile as where β is the dimensionless slope of the pore. If β > 1, then we refer to this as an expanding pore, when β < 1 we have a constricting pore, and when β = 1 we have a cylinder.

Lubrication model for the fluid flow
The pores are typically long and thin and therefore the parameter ϵ 1. For the fluid flow, we solve the lubrication equations, found by considering only the leading-order terms in ϵ in the dimensionless Stokes equations (9), which are given by Equation (14a) yields that the pressure is independent of r, i.e., p = p(z, t). The momentum equation (14b) can hence be integrated in r twice, and after using conditions (2a) and (2b) to fix the constants of integration, we find the following expression for the axial velocity Using the continuity equation (14c), we obtain the solution for the radial velocity, given by We note that the constant of integration in the calculation of u was determined from the symmetry condition (2a). Imposing the no-flux condition (2b) yields that for some function of time F (t) to be determined. Substituting the solution for w in (15) and evaluating the integral in (3) gives the flow rate into the pore, To determine F (t), we assume that there is a constant, known pressure drop along the pore, ∆p = p(0, t) − p(1, t) and use (17) to obtain an expression for the flow rate However, if ∆p is not known, we may alternatively fix the inlet speed along the pore's centreline, w in . In this case, from (15) we find w in = −p z (0, t)h 2 (0, t)/4, which, when combined with (17), gives We will use these results in our time-dependent simulations that follow in the next sections.

Lubrication model for the particle concentration
At leading-order, the dimensionless particle concentration equation and corresponding boundary conditions (10) are given by which can be solved to give that c = c(z, t). This motivates introducing an expansion of the form in which case the leading-order system becomes with boundary conditions Averaging equation (22a) over the cross-section of the pore, i.e., integrating it as follows, allows us to use the boundary conditions (22b)-(22c) for c 1r at r = 0 and r = h, resulting in an equation for the leading-order particle concentration c 0 (z, t) [14]. This needs to be solved in conjunction with the equation for the growth of the pore surface and relevant boundary conditions where the flow rate Q is given at each instant in time by (20).

Numerical simulations
In this section we perform time-dependent simulations of the lubrication model presented in Section 4. The physical parameter values and corresponding dimensionless parameters used in the simulations can be found in Table 1. Figure 2 shows the clogging times (in minutes) as a function of the geometry for Pe = 1, 10 and 1000 (which could be interpreted as variations in the inlet speed or the contaminant particle diffusivity). The mean pore radius is kept fixed at 30 µm, whereas the inlet pore diameter is allowed to range between 3 and 57 µm.  Pore diameter (30-100 µm), R 30 µm Pore length (3-4 mm), L 3 mm Concentration at inlet 100 mol/m 3 Axial speed, W 1 m/s Clogging rate, K 2 3.33×10 −10 m 4 /(s mol) Adsorption rate, For the parameters chosen, the clogging times obtained are consistent with the times observed in experiments, but calibration with experimental data is essential to establish the precise values for the parameters k 1 and k 2 [12]. For Pe = O(1) we find increased clogging times for constricting pores. As Pe increases, better performance is observed for the cylindrical pores. In fact, the clogging times in the large-Pe limit can be estimated from (23b), by arguing that, in this limit, the concentration within the pore does not vary appreciably from the concentration at the inlet. In doing so, we find the clogging time, t * = min(r in , r out )/k 2 , or, in dimensional units This estimate is shown as the dashed curve in figure 2(c). Since good agreement is observed for expanding pores for all Pe, in practice, a pore that widens with depth may be used experimentally to estimate K 2 from (24), given the clogging time and the concentration at the inlet. It is also important to highlight the effect of the parameter K 1 , which is a measure of the rate of particle adsorption at the walls of the pore. We found that as K 1 → 0 the regime of applicability of (24) extends to lower values of Pe. As a result, if the value of K 1 is lowered sufficiently, both the Pe = 1 and 10 cases collapse onto the dashed curve in figure 2(c). This is a consequence of the fact that for lower K 1 the walls interact less with the contaminants, which is also expected to occur when the fluid passes through the pore at high speeds in the high-Pe limit. A measure of the performance of the pore is an estimate of the amount of contaminant removed from the flow until clogging occurs, which can be found from different Pe may not be appropriate in this case. What the calculation highlights, however, is that for lower Pe, better efficiency is achieved for constricting pores than expanding pores, given that for any two configurations that yield identical clogging times, the constricting pore is able to remove a greater amount of the contaminant. This effect diminishes as Pe increases, where the cylindrical pore outperforms the conical configurations.

Non-Newtonian effects
When considering the filtering of rheologically complex fluids, the extra stress tensor is no longer a linear, isotropic function of the components of the velocity gradients, but it requires more complicated constitutive equations. For example, blood exhibits non-Newtonian properties, i.e., shear-thinning, viscoelasticity, thixotropy and yield stress [18]. As a first attempt to study the performance of a pore for non-Newtonian liquids, we consider a simple power-law model based on the Ostwald-de Waele constitutive equation [10], which incorporates shear-thickening and shear-thinning effects. Although its main limitation is that it predicts an unbounded stress in regions where the shear-rate is zero, this model is commonly invoked as a starting point for more complex models, e.g. the Herschel-Bulkley model [18]. The extra stress tensor is defined in terms of the rate-of-strain tensorγ = ∇u + (∇u) T and the shear rateγ = √γ ijγji /2 as where µ(γ) is the (non-constant) fluid viscosity. For power-law fluids, we take where k is the consistency constant (units: Pa s n ) and n the power-law index (dimensionless). In the special case of n = 1 we recover the Newtonian law. For n < 1 we obtain the shear-thinning (or pseudo-plastic) behaviour, i.e., the viscosity is a monotonically decreasing function of the shear-rate. When n > 1 the fluid is shearthickening and the viscosity is a monotonically increasing function of the shear-rate. For such fluids, the equation for conservation of mass remains the same as in (1a), but the momentum equation (1b) becomes with τ given by (25). The boundary conditions, (2), remain unchanged. We non-dimensionalise in the same way as for a Newtonian fluid, as outlined in (7), apart from the pressure, which is scaled by P = kW n L/R n+1 , the corresponding momentum equation under the lubrication approximation becomes supplemented by (14a), (14c) and the same boundary conditions as before. Following the same steps as previously, we find and that (17) generalises to For a constant inlet centreline velocity, the flow rate Q(t) is given by which is combined with (23). Note that when n = 1, (30) reduces to the Newtonian version, (20). Our preliminary studies for non-Newtonian fluids (n = 1) did not yield appreciable differences to the Newtonian (n = 1) case studied in Section 4.3. However, we expect that more complex fluids (e.g. viscoplastic fluids) will exhibit significantly different behaviours, which may be a suitable topic for future investigation.

Filter pores coated with a catalyst
One useful application of the ceramic filters is when the pore surface is coated with a chemical catalyst. The catalyst triggers a chemical reaction involving the contaminant, turning this into a less harmful substance. We assume the phase of the reactant and the product to be gaseous. The relevant equations in this case are those presented in Section 4.2, where now k 2 = 0, since we assume the product does not accumulate on the surface of the pores but is carried away. Furthermore, we study the situation when the system is at steady state. In figure 4, we show plots of the contaminant concentration c versus the depth of the pore z for three different values of pore constriction parameter β defined in (13). In all three plots, we keep the mean radius, R = 1, the same, and we take P e = 0.1, k 1 = 1 for our base case with cylindrical pores.
The best removal efficiency is achieved by the constricting pore (β < 1), identified by the lowest concentration of contaminant at the outlet. This is due to the fact that a constricting pore has a higher surface area near the top, where the contaminant concentration is highest. As the contaminant concentration falls due to its removal from the flow, the pore area available for chemical reaction reduces accordingly. We also see that for the constricting pore the removal rate is more uniform (i.e., the drop in contaminant is closest to linear) than the other two configurations.

The effect of number of pores
All of the above analysis presented so far is for a single pore. To model a membrane, which consists of multiple pores, it is important to explore the effect of the geometry on the possible distribution of these pores. To this end, we consider a very simple illustrative comparison between a 2-cylinder-pore (β = 1 in (13)) and a 3-cone-pore (β = 0.9) configuration (see figure 5). We keep the same cross-sectional area of the filter and total influx of contaminant, Q, in both cases. Thus, we have where Q 2 and Q 3 are the corresponding fluxes through the pores in the 2-pore and 3-pore configuration, respectively. In figure 6, we show plots of the contaminant concentration versus the depth of the pore for the two configurations. The difference between the performance in the two cases is minimal, with a slight advantage for the 3-pore case over the 2-pore case. Further investigation is necessary when arbitrary numbers and shapes of pores are used. In addition, any structural limitations of the filter, i.e., how close the pores can be to each other, must also be taken into account.

Introduction
So far we have performed simulations that exploit the slenderness of the pore geometry by employing a lubrication approach. In this section we relax the assumption of pore slenderness, and use two-dimensional simulations to model the changes in the pore geometry as contaminant particles are absorbed onto the pore walls. We investigate the shape of the pores that provides the longest lifetime of the membrane before the pore clogging occurs. In this section we consider two-dimensional Cartesian configurations rather than axisymmetric.
We consider the dimensionless model discussed in Section 3. Since we do not make any assumptions on the pore size or shape, we simply choose the lengthscales as L = R = (r in + r out )/2 and the velocity scales as U = W = w in (0), where w in (r) is the inflow velocity profile. Then, we supplement the fluid-flow problem (9) with the following boundary conditions: where n is the outward unit normal to the boundary of the computational domain, ∂Ω (see figure 1). Equation (32a) is a Dirichlet boundary condition for the inflow velocity, while (32b) enforces no stress on the outflow boundary. Similarly, equation (10) for the concentration of contaminant is supported by the following initial condition and boundary conditions on the inflow and outflow boundaries: The initial condition for equation (11), which determines the shape of the pore wall before contaminant adsorption occurs is given by (13).

Free-boundary problem discretisation
The final system of equations consists of (2), (9)-(11), (32)-(33) and represents a free-boundary-value problem: the contaminant, by depositing on the free surface h(z) as time progresses, modifies the geometry of the domain considered. This would force a re-discretisation of the domain at every time step. To circumvent this issue, we treat modifications to the domain geometry implicitly, by considering a simple linear map F from a simplified reference domainΩ, (a unit square), to the deformed physical computational domain Ω. This map, whose action is sketched in figure 7, is defined as follows: , with Jacobian J F = ) .
Here and below, we mark with a hat all variables and operators defined in the reference domain.
Introducing F allows us to solve the target equations on the fixed reference domainΩ, encapsulating modifications to the geometry directly within the map. In order to take into account the effect of these modifications, though, the target equations need to be modified accordingly. The variational form of the flow problem in the reference domainΩ reads ∫Ω (∇û whereŝ is the test function. We discretise the problem using finite-element library FEniCS [1]. A mixed finiteelement method is employed to solve the flow problem with linear and quadratic Lagrange elements for the pressure and velocity, respectively. The mass-transport problem is discretised using quadratic Lagrange elements.

Pore clogging
The purpose of the numerical experiments conducted in this section is to investigate how the clogging occurs for pores with different initial shape and under different operational conditions. Table 2 shows the set of dimensionless parameters used for the simulations.   We find that the pore shape with r in = 1.25 and r out = 0.75 provides the best usage of the pore space and, when the pore is clogged, the whole pore volume has been used for the contaminant adsorption. On the other hand, the other two pore shapes have a lot of unused pore space left after the clogging has taken place. This arises as a result of the fact that the contaminant concentration decreases with depth due to adsorption. The final times also indicate the preferential shapes: t = 1.72 for the pore with the larger inflow radius than outflow, t = 1.36 for the straight pore, and t = 1.02 for the pore with the smaller inflow radius than outflow.
A faster flow regime, Péclet number Pe = 10 3 , is considered in figure 9. Similarly to the previous case, we show the concentration distributions with the velocity streamlines for the three different shapes at initial and final times. In this case, we observe that the best performing pore shape is the straight one with r in = r out = 1, which has the final time t = 1.26, while the other two pore shapes both last until t = 0.96. This is because at higher flow rates, the concentration of contaminant is approximately constant throughout the pore and does not fall with depth as was observed in the slower-flow case. Thus we conclude that the optimal shape depends on the operating conditions.
As future work, one should consider other performance characteristics as well, such as pressure drop and the efficiency, in addition to the pore lifetime before clogging.

Effect of flow angle
In this section, we investigate how the angle of the fluid velocity at the inflow affects the pressure drop and the membrane performance with different pore shapes. Here we consider only the initial membrane performance, namely before the shape of the pore starts changing due to the contaminant adsorption. Then, the problem consists of equations (1), (2b), (4), (5b) and dimensional form of the inflow and outflow boundary conditions (32) and (33). To perform the simulation, we use an in-house software tool developed and discussed in [13].  Inflow velocity magnitude |u| 10 cm/s Inflow concentration c in 10 14 particle count/m 3 Adsorption rate K 1 10 −5 m/s Viscosity µ 1.85 × 10 −5 Pa s Porosity ϕ 60 % Pore height L 200 µm Inlet radius r in 40, 60, 80 µm Outlet radius r out 80, 60, 40 µm Table 3: Input parameters for 2D simulations of different inflow angles.
We consider seven different inflow angles α = 0°, 10°. . . 60°. Figure 11 shows the velocity streamlines coloured with the velocity magnitude with the inflow angle α = 30°for three different pore geometries, which we will refer to as constricting, straight or expanding. Corresponding concentration profiles are shown in figure 12 while figure 13 shows the pressure drop. We define the initial efficiency via and plot this in figure 14 for different inflow angles for the three types of the geometries. We observe that the straight pore has the lowest pressure drop, while the expanding pore has the highest one. On the other hand, while having the intermediate pressure drop, the constricting pore provides the largest filtration efficiency. The straight pore exhibits the second best efficiency. The worst performance is displayed by the expanding pore in terms of both the pressure drop and the filtration efficiency. For the different inflow angles, we observe that pressure drop decreases and the efficiency increases slightly as the inflow angle increases. In summary, the numerical experiments show that the straight and constricting pores exhibit promising behaviour in terms of the pressure drop and the efficiency especially when applying inclined inflow velocity. In future work, using numerical simulations one can further investigate how the inclined velocity influences particles deposition on top of the membrane and on the pore walls and how it can be used to clean the contaminant fouling on the surface of the membrane. Furthermore, to obtain a more accurate picture of this kind of filtration experiments the problem needs to be extended and considered in 3D.

Discrete model
All of the models we have analysed so far (mostly) consider a single pore, and assume that the overall filter behaviour is simply given by multiple pores all behaving in an identical manner. For small particles, this assumption is valid, but for larger particles, where only a small number of particles are needed to block a pore, this breaks down. In figure 15 we show a schematic diagram of the situation we consider.  In our model, we assume that the simple membrane has uniform thickness, which comprises an array of m regularly spaced uniform holes of radii r 0 , called pores. These pores run from one side of the membrane to the other. This type of structure, for example, describes track-etched membranes, where the pores are cylindrical or cone-shaped (see for instance [2]). As the fluid passes through the pores of the membrane the particles are potentially retained with some probability. The retention mechanisms are rather complex and depend on properties of the membrane and the particles themselves, and they could even change during the filtration process. In what we consider here, we assume that the probability of a particle sticking within the membrane is 1, so every particle is retained by the membrane.  A stochastic simulation is proposed in Griffiths et al. [9]. The authors consider a three-step blocking process in a constant-pressure filtration setting. In each pore, uniformly sized particles first constrict the pore by adhesion on the pore wall. Then, as the radius of the pore decreases, they partially cover the pore and finally form a layer of particles on the membrane surface, termed a filtercake. The same blocking process occurs in all of the pores, and therefore, the flux through a given pore is only dependent on the number of particles it has retained. However, crucially, each pore may contain a different number of particles. In [9], multiple stochastic simulations were run to compute the flux decline for the entire membrane. The probability of a particle landing in a pore at a given iteration step depends on the flux through that pore, and the flux through the entire membrane is obtained by summing over the fluxes of the individual pores.
The number and order of arrival of the retained particles for a given pore represents the state of the pore. In this section we adopt an alternative approach and calculate the expected flux through the membrane. Our analytic model based on expectation values extracts the correct smooth result that would be observed experimentally.
We begin by forming the so-called state matrices, denoted by g, which represent the state of the system, i.e., the pore in which each particle goes in, and the order in which this happens. These state matrices have binary entries. Each row of the matrix g represents a particle and each column of matrix g represents a pore. We assume that there are m pores and n particles, so there are a total of m n distinct state matrices. For instance, for the case m = 3, n = 4, we can have a state matrix of the form g =     g 11 g 12 g 13 g 21 g 22 g 23 g 31 g 32 g 33 g 41 g 42 g 43 which implies that the first particle goes into pore 1, the second particle goes into pore 2, the third particle goes into pore 3 and finally, the fourth particle goes into pore 1 again. From the nature of g, we can infer that there can be at most one 1 in each row. Next, we form the flux matrix, M . Each row of M represents the flux after the jth particle has entered a pore, which we define as Q(j). For a system of three pores, the flux matrix is then of the following form . . .
and so, the expectation of flux then after l particles is given by where subscripts denote the entry of that vector. We claim that it is sufficient to consider only a few pores to provide a good prediction and so we only consider three We consider a flux that has the following form Q = [10,9,7,3,0]. The red line depicts the flux as we would naively expect it to decline if the filter blocked by first filling each pore with a single particle before depositing a second particle in each pore, and so on, i.e., Q exp = [30, 27, 21, 9] for n = [0, 3,6,9]. The blue line shows our expectation-value prediction using the procedure described above, which already expresses a smoothed-out version despite us only considering three pores.
pores for our simulations here. In figure 16 we consider a scenario where the flux in a pore is defined by Q = [10,9,7,3,0], so that the flux drops by one unit when the first particle is accepted, a further two units when the second particle is accepted, a further four units when the third particle is accepted, and a further three units when the last particle is accepted. We present here the naive prediction if we assumed that the particles were to deposit in the pores in a uniform manner, i.e., one particle occupies each individual pore before the pores accept a second particle, and so on. This prediction (shown in red in figure 16) shows clear jumps in the predicted flux decline, which would not be observed in practice. Our expectation value approach, despite only here applying to a membrane comprising three pores, already predicts how the real experimental curve would appear smoothed out. While we have only presented a simple model for the expectation value of the membrane at a given instant in time, this model demonstrates clear potential in the ability to predict the smoothed flux versus throughput curve that would be observed experimentally. Increasing the number of pores considered would improve the smoothness of this prediction.

Conclusions and recommendations
In this report we have considered a series of distinct yet complementary models to gain insight into the behaviour of the SSL ceramic membranes. Here, we summarise our conclusions: • By exploiting the slenderness of the pores in the membrane we may make a lubrication approximation to the governing equations. The resulting model is amenable to analytical techniques, and allows us to predict how the flux through the membrane will vary with pore profile, (19).
• Numerical simulations in the lubrication regime provide a prediction for how a pore will clog due to contaminants adhering to the pore walls. We can predict the time to clogging of the membrane as a function of the pore angle. We find that a membrane comprising pores that constrict with depth takes longer to clog (figure 2).
• Numerical simulations in the lubrication regime provide a prediction for how much contaminant is removed before the membrane clogs. We find that a membrane comprising pores that constrict with depth removes the most material ( figure 3).
• We can generalise the lubrication model to add non-Newtonian effects. While we presented the framework here, more work is needed to study the behaviour in a particular case.
• Numerical simulations in the lubrication regime provide a prediction for how the filter behaves when removing contaminant through a catalytic reaction (which does not clog the membrane). We find that a membrane comprising pores that constrict with depth removes the contaminant in the most uniform manner with depth ( figure 4).
• Numerical simulations in the lubrication regime provides insight into how the arrangement of pores affects the removal efficiency (figure 6). A more comprehensive parameter sweep is required to understand the role of pore arrangement fully.
• We may relax our assumption of slender pores and derive a fully 2D framework.
To account for the changing pore shape due to blocking we use a movingboundary mesh formulation. The results provide detailed insight into how the flow behaves in the pore and how the pore evolves due to surface deposition at both low speeds (figure 8) and high speeds ( figure 9).
• The full numerical simulation also allows us to explore the effect of flowing the contaminant feed at an angle to the membrane surface. We find that the pressure drop is highest in a constricting pore and lowest in a straight pore, irrespective of the inflow angle. A widening pore has an intermediate pressure drop ( figure 13). The efficiency is also highest for a constricting pore. However, the widening pore has the lowest efficiency and the straight pore has an intermediate efficiency ( figure 14).
• Finally, we presented a probabilistic model based on expectation values for the membrane as particles deposit into the filter. The model is able to provide analytic predictions for how the flux through the filter falls with time, capturing the variability in the number of particles that are contained within each pore in the membrane, without the need for time-consuming stochastic Monte Carlo simulations ( figure 16).
In terms of next steps, we have: • Developed a 10-week mini-project in the Industrially Focused Mathematical Modelling (InFoMM) Centre for Doctoral Training (CDT) at the Mathematical Institute, University of Oxford which followed on with this work.
• Submitted a proposal and secured an EPSRC Impact Acceleration Account Award at the Mathematical Institute, University of Oxford, to explore in more detail the probabilistic and numerical models presented in this report.
• Had preliminary discussions on modelling aspects of this challenge with non-Newtonian flows, at the School of Mathematics, Cardiff University.