Optimizing the Axial Resistance Proﬁle of Submerged Hollow Fiber Membranes

: Submerged hollow ﬁber membranes (SHFMs) are used for a wide variety of applications. Their applicability and their advantages, however, depend crucially on the prevailing hydrodynamics within single ﬁbers. In this respect, the non-uniform distribution of transmembrane ﬂux is a known problem related to inhomogeneous membrane fouling and disadvantages for cleaning. To address this problem, we propose an approach to homogenize transmembrane ﬂux by varying the local membrane resistance using optimal control methods for the ﬁrst time in SHFM research. Based on an established model, different scenarios are optimized, namely with different ﬁber lengths and inner radii. In addition, a double-end setup is explored. It is shown that the optimization goal is reached very well in all tested cases, which underlines the general validity of our strategy. Further uses and extensions of the optimization method are provided, as well as hints for the practical implementation of the suggested measures.


Introduction
Submerged hollow fiber membranes (SHFMs) are a versatile separation technology. SHFMs are employed for various filtration problems such as those encountered in water treatment, e.g., microfiltration, ultrafiltration, and reverse osmosis, or in membrane bioreactors. Compared to flat sheet or tubular membrane modules, SHFMs have the advantages of high packing density, good cleanability, small energy expenditure during operation, and low fabrication costs [1][2][3][4]. SHFMs consist of bundles of fibers, usually called modules, which are either horizontally or vertically submerged in a pressurized liquid reservoir. If the fibers are not submerged in a tank, but incorporated in a cartridge with a cylindrical casing, the membranes are commonly only referred to as hollow fiber membranes; however, the prevailing physics remain largely the same. In most setups, a suction pump is used, often in the mode of constant flow, to transport the permeate from the membrane module [2]. Importantly, the said advantages, mainly concerning cleanability and energy expenditure, crucially depend on the prevailing hydrodynamics. In this respect, the arrangement of single fibers within a module has been found to be important [4,5]. However, particular emphasis has also been placed on the hydrodynamics in single fibers. The non-uniform profile of transmembrane flux, with the highest flux at the membrane outlet, has been repeatedly claimed as a significant disadvantage [1][2][3][4]. When operated in non-cross-flow mode, the profile of transmembrane flux leads to a locally varying fouling behavior. In the cross-flow mode, where often bubbles are used to enhance cleaning in case of vertically oriented SHFM modules, it is reported that below a critical flux, no problems occur because the influence of cross-flow is pronounced enough to prevent significant clogging. However, due to the higher values of transmembrane flux at the outlet, this critical flux is easily superseded there [1,2,6]. In previous work, the influence of basic design parameters, such as fiber length and fiber diameter, on the profile of transmembrane flux was studied [1][2][3].
Our aim is also to address the known problem of a non-uniform profile of transmembrane flux. However, the problem is addressed differently than before: we propose to use another influencing factor as the control variable, namely the local membrane resistance, which is varied along the fiber's length coordinate. Except in the conclusion, the practical implementation and questions of manufacturing are not considered in this article; rather, the focus lies on what one would ideally like to manufacture. The introduced method is a model-based optimization approach using optimal control theory. For the model, we rely on the well established formulation by Chang and Fane [1]. The same model is also used in later publications of their group where it is employed for parametric studies [2], extended to also include clogging [7], and applied to aid the interpretation of experimental findings [8]. The same model is also utilized in various other works [9,10]. Based on a slightly different derivation, Ye et al. [11] arrive at a very similar mathematical formulation, which also holds true for the model of Liang et al. [12]. In newer modeling approaches, e.g., by Li et al. [3], the model by Chang et al. is still used as a reference case and very similar results are achieved compared with the new model. For this reason, the model by Chang and Fane [1] is considered a well-tested standard case and a there is a good basis for incorporating it into our optimization approach.
As a reliable mathematical process description exists, the focus of this work is not on modeling but on proposing a new optimization method, which is a common goal in process systems engineering [13]. Even though different works have aimed at optimizing SHFMs [2,[14][15][16][17], this is the first study where optimal control theory is applied to hollow fiber membranes. In general, dynamic programming or optimal control methods are only rarely encountered in filtration research. Previous investigations in the context of membrane science were done by Blankert et al. [18,19], Zondervan et al. [20], and Paulen et al. [21,22], who determined the optimal time-trajectories for operation; in addition, some newer work addresses time-optimal control in the field of filtration [23]. In contrast, optimal control is conducted here along a spatial dimension, an aspect only scarcely addressed within classical optimal control research [24][25][26]. Therefore, this can be considered an innovative methodological twist, which also implies that optimal control is applied for the design of hardware instead of, as usual, for optimization of the operation of given equipment. The presented approach builds on previous work by the authors where depth filter media were optimized along their flow direction, i.e., also along a spatial dimension [27][28][29]. Therefore, the novelty claim of this work does not consist of a new optimal control algorithm, but rather in its somewhat unusual utilization by optimizing the control variable in a spatial dimension and mainly in employing optimal control for the first time in SHFM research. Thus, when we speak of a "new optimization method", a new tool for optimizing the design of SHFMs is meant.

Theoretical Background
As discussed, we will base our work on the established model presented by Chang and Fane [1]. They, in turn, rely on earlier work, mainly on Apelblat et al. [30], who derived a model for fluid exchange between blood capillaries and the surrounding tissue through permeable membrane walls. Doshi et al. [14] and Bruining [31] already used similar equations to predict the hydraulic behavior of hollow fiber membranes. These modeling approaches consider only the hydraulics of single fibers as also illustrated in Figure 1, where a basic model sketch is provided. In the approach by Chang and Fane [1], two main effects are combined: laminar pipe flow and Darcian transmembrane flux. The model derivation is briefly recapitulated to prepare the model for its use in our optimal control approach. Volumetric flow rate Q in the cylindrical hollow membrane is described by the Hagen-Poiseuille equation: where r is the inner radius of the membrane, µ the fluid's dynamic viscosity, p the pressure inside the membrane, and x the spatial coordinate, starting at the dead-end side of the membrane. Transmembrane flux J is modeled based on Darcy's law: ∆P is the local pressure difference across the membrane's surface, R is the local membrane resistance, and p o is the fluid pressure outside of the membrane. Please note that in normal operation R is the sum of the resistance of the clean membrane R m and some additional resistance R d due to deposition of impurities on the membrane surface, i.e., R = R m + R d ; this will be addressed later on in this section. The local mass balance describes the increase in flow rate Q along the length coordinate x due to transmembrane flux J. If Equations (1) and (2) are introduced into this mass balance, one obtains where λ is defined as to abbreviate further equations. Considering only the pure hydraulics of submerged hollow fiber membranes, i.e., without any effects of deposition or fouling on resistance (R = R m ), Chang and Fane [1] derived an analytical solution for Equation (4). Assuming no axial inflow, i.e., dp/dx = 0 at x = 0, plus providing the average transmembrane flux J, i.e., the resulting integration constants can be determined and the solution becomes [1]: or, using Equation (2), where L is the length of a single fiber. However, three things are important to note with respect to these results. First, inclusion of hydrostatic pressure for vertically oriented fibers, as done by Chang et al. [1,2], does not influence the overall behavior. Hydrostatics affects the fluid surrounding the membrane, i.e., p o , as well as the fluid column inside the membrane, i.e., p. Therefore, it does not impact the flux J(x). Equation (8), therefore, has the same form for vertically and horizontally oriented membrane fibers. However, hydrostatic pressure can be easily included if the absolute pressure values inside the membrane or on the outlet of single fibers is of interest. Secondly, the classical model describes steady-state behavior. Only if fouling is included can the behavior become transient due to the fact that more deposits are accumulated over time on the membrane fiber. As our optimization approach only addresses the hydrodynamics of single fibers without fouling, we are also dealing with a steady-state problem. The third point we want to stress is the limited validity of Equations (7) and (8). For obtaining these analytical solutions, it is implicitly assumed that the resistance R does not change with x, as can be seen from the foregoing derivation. This is, however, no longer true when fouling is considered because the non-constant transmembrane flux J(x) will also cause a locally varying deposition of impurities. We assume that a constant R is also no longer feasible when the membrane resistance R m is made the control variable in our optimal control approach because it is then, by definition, dependent on the spatial variable x, i.e., R m (x). For these reasons, Equations (7) and (8) cannot be used for the main part of this work and we need to rely instead on solving Equation (4) directly by numerical methods. As already explained, when the deposition of impurities on the membrane is considered, resistance R is modeled as R m + R d , i.e., as the sum of the pure membrane resistance R m and the additional resistance due to deposition R d . R m is a property of the membrane and, therefore, independent of time. However, it may be dependent on space, i.e., R m (x) as in our optimization. Optimization of R m (x) is thus a typical design problem. In contrast, R d is a function of space and time, i.e., R d (t, x), as it is not a function of design but of operation. For reasons of readability, we will omit these independent variables in the following, and include them only if special attention is drawn to them. R d depends on the amount of liquid w that passes through the membrane per unit filter surface: which influences the local resistance R d in the following way: where C and n are constants describing the increase in resistance for different filtration mechanisms. Equation (10) was introduced by Blankert et al. [18] and is intended as a modern formulation of the classical laws of filtration proposed by Hermans and Bredée [32] and Hermia [33]. In these laws of filtration, the exponent n indicates different filtration mechanisms; it is 0 for incompressible cake filtration, 1 for a mechanism referred to as intermediate blocking, 3/2 for standard blocking, and 2 for complete blocking [18,34]. Deposition on hollow fiber membranes has been already investigated in the literature [7,9,35]. Lee et al. [35] found that cake filtration and standard blocking are the dominant fouling mechanisms for hydrophilic and hydrophobic membranes, respectively. Liu et al. [9] included cake filtration and also modeled cake compression. Please note that Equations (9) and (10) are only unrestrictedly valid for a non-cross-flow operation mode of SHFMs. Membrane fouling is only briefly treated here in order to motivate our optimization approach and aid interpretation of the results. We will only consider cake filtration, i.e., n = 0. In this case, C is related to the resistance of the filter cake formed by deposited impurities [18]. Please note that for pronounced cake formation, the one-dimensional Cartesian description that has been used can also be easily transferred to cylindrical geometries; various formulations are found in the literature [36][37][38]. Having mentioned this, we will not further consider the cylindrical fiber geometry and treat the deposition as a pure one-dimensional Cartesian problem as it is also approached in the cited literature on modeling of hollow fiber membranes; this is a valid approximation for small cake heights. Even though only cake filtration is treated, our reasoning and the optimization results, however, are valid for all other filtration modes as well. This is due to the fact that for all filtration mechanisms, the local fouling state depends on w(t, x) which, according to Equation (9), is the time integral of J. Thus, if our optimization goal of a constant local flux J(x) is reached, single membrane fibers equally exhibit constant local values of w(x) and, therefore, a constant resistance R d (t, x) for each point in time t. For this reason, cake filtration serves only as one example of fouling. The reasoning, however, remains general.
With respect to the overall approach, please note that more complicated models have also been proposed in the literature. One obvious disadvantage of the model based on Equations (1)-(3) is that it relies on the simplest coupling of transmembrane flux and pipe flow. Physically, however, the ideal axial pipe flow may be influenced by the perpendicular flow through the membrane wall. For this reason Kim et el. [39] and Li et al. [3] proposed alternative models that take this effect into account. Nevertheless, we will use the basic modeling approach of Chang and Fane [1] as this model has proved suitable for many applications, especially when low transmembrane fluxes are encountered, and can be still solved with reasonable computational efficiency required when addressing optimal control problems.
All parameters needed to solve the introduced model equations are summarized in Table 1. Single parameter values are varied in different scenarios considered in the results section; these values are shown in the last row of the table. The values for the reference case are taken from Chang et al. [1,2] except for the value of C, which is chosen so that typical time scale and fouling effects from the literature are met [1]. In addition, other publications show that the parameter values used in this study are in the typical range [3,8,11,31].

Numerical Methods
All numerical simulations and optimization operations are performed using MATLAB (version: 2019a, supplier: The MathWorks, Natick, MA, USA). As explained, the basis for solving the model is Equation (4). To meet our reference case of Chang et al. [1,2], however, the simple solution of this second-order differential equation must be supplied with an additional building block. This is due to the fact that for solving Equation (4), boundary conditions for p and dp/dx are required. Chang et al. as well as other groups [1,2,35] use the average transmembrane flux J, as defined in Equation (6), as an overall constraint. However, when the analytical solutions of Equations (7) and (8) are no longer valid, e.g., due to locally varying resistances, J can only be computed from the final simulation results. For this reason, a self-implemented bisection method is used to determine the boundary condition for p at x = 0 so that the desired value for J is met. Bisection is conducted until the difference between upper and lower interval becomes smaller than 10 −3 times the mean value of the solution. This threshold proved sufficient to achieve inaccuracies of less than 0.02% in all performed simulations, i.e., J deviated less than this percentage from its predefined value. The procedure is insensitive to the chosen start values. A validation of this procedure is provided at the beginning of the results section. The said bisection method is used iteratively with the MATLAB solver ode45, by which the two coupled ordinary differential equations are solved; ode45 is used with its default settings. The boundary condition for dp/dx at x = 0 is 0 when no in-or outflow at the lower end takes place. The situation changes in a double-end setup, as explained later; see Equation (12).
Membrane fouling due to cake formation is computed by an explicit method, i.e., integration moves forward in time without any inner iterations. The local flow profile is obtained by the approach just described. Current local resistance R(t, x) is determined according to Equation (10) based on the local cumulative volume w(t, x), as defined in Equation (9), for each point in time. Initially, the membrane is assumed to be unclogged, i.e., R d (t, x) = 0 at t = 0. Using the profile of R(t, x) at time T, J(t, x) for T + ∆t is computed and so on. A time-step size ∆t of 30 s is used, which is chosen according to a convergence study, i.e., a finer time discretization does not change the results anymore.
As mentioned, we aim at a constant transmembrane flux. This leads to the following cost functional used in the optimal control computation i.e., the variance of the flux; please note that Equation (11) is a Lagrangian-type objective functional. Thus, the optimization task is min( f ), which is achieved by optimizing the path of R m (x). Note that this is a time-invariant problem as we only optimize the membrane resistance R m as a function of the fiber's length coordinate x. Fouling is not included in the optimization procedure. However, we argued on physical grounds in Section 2 that a constant local flux also leads to a constant deposition along the fiber length. Numerically, optimal control is approached by a direct single shooting algorithm based on discretization of the control variable [40][41][42], i.e., R m in our case. This is a so-called direct solution method where the cost functional is minimized straightforwardly, in contrast to indirect methods where a Hamiltonian together with necessary and, if possible, sufficient conditions is employed to determine the optimal paths [41]. Initially, the unknown trajectory of R m (x) is approximated by a linear profile that is optimized. The optimal linear profile, defined by two points, is intersected by a third point in the next iteration, for which the previous linear solution is used as the initialization. In all further optimization steps, additional points are introduced between previously optimized locations; all points are linearly connected. This iterative division and optimization procedure is described in more detail in other publications by the authors [27,28]. Numerical optimization is conducted using MATLAB's function fminsearch with its default settings. The sought-for optimal path is approximated by six iterations, corresponding to 33 points, which yields sufficiently smooth profiles in the current case, as can be seen in the results section.

Results and Discussion
First of all, a validation of our numerical solution method, described in the last section, is provided. For this reason, the analytical solution of Equation (8) is compared with the numerical solution of Equation (4). As can be seen in Figure 2, both results show perfect agreement for the parameters of the reference case, as given in Table 1, and which also remains true when these parameters are varied. Therefore, our solution procedure, including the described bisection method, is seen as valid. Due to the locally varying transmembrane flux, a locally varying fouling behavior can be expected in non-cross-flow mode. This behavior was also found in the literature [43] and is, as described, seen as a problem in the operation of submerged hollow fiber membranes. Figure 3 shows the increase in local resistance R(t, x) and its transient influence on local transmembrane flux J(t, x). Increase in local resistance is most pronounced at the membrane outlet (x = L) where the transmembrane flux is highest, as already shown in Figure 3a. The coupling of fouling and flow is due to the fact that each fluid element also contains some amount of foulants and that, therefore, a higher amount of liquid than passed locally through the membrane leads to a more pronounced clogging at this location. Fouling, in turn, causes a local decrease in flux as depicted in Figure 3b. Therefore, a negative feedback loop is created, i.e., at the fiber outlet, flux is highest, which is also related to the most pronounced clogging which, in turn, causes the strongest decline in local flux. Nevertheless, deposition of impurities on the membrane surface is still very inhomogeneous as can be seen from the changing resistance profile over time in Figure 3a. As explained in the introduction, inhomogeneous clogging behavior is disadvantageous for membrane cleaning and it is, therefore, desirable to avoid it. Due to the coupling of fouling and flux, a constant transmembrane flux also implies constant fouling along the membrane length, as argued in Section 2. Figure 3, therefore, retrospectively justifies our chosen optimization goal. The said behavior is especially important in the non-cross-flow mode, but when cross-flow is applied it is also required that the critical flux is nowhere superseded, as explained in the introduction. Therefore, a con-stant flux with a suitably chosen value is the safest way to assure that the critical flux is not reached.
In our first optimal-control case study, our reference case by Chang and Fane [1] is optimized to achieve a constant flux. Please note that here and in all of the following optimization scenarios, we deal again with problems of membrane design that are therefore time-invariant. All shown resistances are membrane resistances R m . In Figure 4, the optimized solution is compared to the gradient-free membrane with the same average flux J. It can be observed that the optimization goal was ideally reached, meaning that the local transmembrane flux is indeed constant along the membrane length. This is a result of the optimized profile of local transmembrane resistance R m (x), which is lowest at the inlet (x = 0) and highest at the membrane outlet (x = L). In order to investigate the influence of different membrane parameters and to underline the general validity of the proposed method, two scenarios with geometrical variations are considered. First, a membrane with an increased length is addressed; second, a membrane with a larger inner radius is studied. These parameters were also varied in the literature with the goal of investigating their influence on the axial profile of transmembrane flux [2]. Figure 5 shows the corresponding results for the longer membrane fiber (L = 1.5 m), both optimized and gradient-free. For achieving the same overall flow rate per fiber as in the reference scenario, which is obtained by multiplying the average flux with the total surface of the membrane fiber, an average flux of J = 20 L m −2 h −1 is used here. The findings are juxtaposed with our reference case of a one-meter-long fiber. It can be seen that longer membranes still cause more inhomogeneous flux profiles, which need to be compensated by a more pronounced variation of local membrane resistance. However, using the optimized resistance profile also given in Figure 5, a constant flux can be achieved. Increase of the inner membrane radius has the opposite effect compared to increasing the length of fibers, i.e., it homogenizes the flux profile. In Figure 6, results for a fiber of 1 m length but an inner radius of 0.3 mm are shown, i.e., 50% larger compared to the reference case. In order to keep the total flow rate constant, again an average flux of J = 20 L m −2 h −1 is used. The gradient-free membrane is compared to the optimized gradient membrane and it can be seen again that the optimized membrane has a completely homogeneous transmembrane flux. It could be argued that membranes with such strong deviations in local resistance might be unrealistic, a point we will also address in the conclusions section. In response to this issue, we also propose a method by which a homogeneous transmembrane flux can be achieved without such strong gradients of R m (x). To achieve this, we turn to a double-end scenario where suction is applied at both ends of the membrane, as also encountered in the literature [4,10]. This setup leads to the following boundary condition at x = 0: dp dx = − 8µQ where s is the share of the total flow rate that is extracted at x = 0. Positive values of s correspond to inflow at x = 0, negative values to outflow. For the previous results, s = 0 was used as no inflow or outflow being assumed at this location. Now, s is set to −0.5, i.e., half of the overall flow rate is extracted at x = 0. In Figure 7, the optimized resistance profile as well as the corresponding flow profiles are shown; both are compared again to the gradient-free scenario. Please note that for these results forward computation was performed directly; however, for determining the optimal control solution, the domain was split, i.e., left and right sides of the curve were computed separately for membrane segments of 0.5 m each. This was done to circumvent possible numerical problems. When splitting the domain, s is again set to 0 for one fiber end (as in all cases above). Both halffiber solutions are subsequently stitched together at the s = 0 location. It can be seen from the solution that the chosen strategy is valid, as the overall optimization goal is perfectly reached and the prescribed flux is met as well. Additionally, Figure 7 illustrates that the double-end setup indeed requires only the smallest variation of membrane resistance compared to all previously studied cases in order to achieve a constant transmembrane flux. This can be explained by the fact that even without optimization, the double-end mode already has a more homogeneous flow profile compared to the reference case. For further comparison of the four different scenarios, each for a homogeneous, i.e., gradient-free, membrane and an optimized gradient membrane, the inhomogeneity of the flux profile is evaluated; Table 2 summarizes the outcomes. As a measure for inhomogeneity of flux we rely on f / J 2 , i.e., the variance of flux f , as defined in Equation (11) scaled by the square of the average flux. The chosen scaling makes sense because in the nonoptimized cases f is in the order of J 2 . As already observed earlier in the article and in the literature [2,4], an increasing length of the membrane fiber leads to a larger inhomogeneity of the flux profile whereas a larger inner radius decreases flux inhomogeneity. In addition, a double-end setup with the same throughput has a lower inhomogeneity of flux. Table 2 shows again in a quantitative manner that the optimization goal is reached very well in all investigated cases, i.e., that all optimized membranes have an ideally homogeneous flux profile. Relying on the previous results, it is easy to show that all optimized membrane fibers also show a completely homogeneous fouling behavior in the non-cross-flow mode. In contrast to the outcomes for the reference case presented in Figure 3, where a pronounced difference in local deposition is visible, a constant transmembrane flux also causes a constant deposition profile on the membrane surface. Note that this is true not only for the mechanism of cake filtration, but also for all other filtration mechanisms discussed earlier because w(t, x), as defined in Equation (9), does not differ along the fiber length. It is expected that this has a positive influence on membrane cleaning. In addition, when cross-flow is applied, the homogeneous flux profile in optimized membrane fibers ensures that the critical flux is nowhere superseded.

Conclusions and Outlook
In this article, we introduced a new optimization method to determine profiles of local resistance which lead to a constant transmembrane flux along the length of submerged hollow fiber membranes. The method is based on optimal control theory and relies on the established and repeatedly validated model by Chang and Fane [1]. Using different scenarios with varying fiber lengths and radii, the broad validity of the approach was illustrated. In all cases, the goal of a locally invariable flux was perfectly reached. In addition, a double-end setup was investigated. This setup already has, per se, a lower inhomogeneity of transmembrane flux, which implies that only a smaller variation of local resistance is required in order to achieve a completely homogeneous profile of transmembrane flux.
As it might be objected that locally varying the resistance of submerged hollow fiber membranes is hardly feasible from a practical point of view, the double-end setup might be a good starting point for an implementation in practice. However, in general, we see a change in filtration research from what we call the "microtheoretic paradigm", where the focus of research turned to the theoretical understanding of pore-scale effects in filters, to a research phase that can be termed the "micromanipulative paradigm" [44]. Now, the microscale of filters can be influenced in previously unimaginable ways. A key role is played here by different 3D-printing and additive manufacturing technologies [45][46][47][48]. In addition, local and selective thermal or chemical treatment on very small length scales can be mentioned. Based on this background, we can also think of ways to locally influence the resistance of single membrane fibers. One can imagine, e.g., a locally varying coating or a targeted deposition of micro-or nanoparticles on the membrane surface, which could be fixed by sintering processes. Influencing the local pore size, which directly impacts resistance, by thermal or chemical treatment might also be an option. A variation of the local lateral thickness of the fibers, as already applied in ceramic gradient membranes [49], might be another means to achieve varying resistances. Obviously, it is also possible to combine some of these manufacturing processes to realize the required gradients.
Besides these direct practical uses of the obtained results, the method itself might be a beneficial basis for further use and development. With the same basic approach, other optimization goals can be addressed, e.g., the resistance profile can be determined so that only a certain maximal flux is not superseded, meaningfully the so-called critical flux. Whereas our algorithmic take on optimizing the axial resistance profile proved sufficient for the goal of a constant flux, other optimization goals might result in more complicated problems with their own challenges concerning numerical stability. In this respect, a wavelet-based optimal control could be promising [50,51]. Besides its use for single fibers, our method can also be combined with design methods for whole modules. These methods account for the influence of fiber position on the profile of transmembrane flux, which also influences the optimization goal for single membrane fibers. Such a problem could be addressed by combining computational fluid dynamics (CFD) simulations of whole modules with our optimization method for single fibers. Additionally, the method might be adapted to other membrane types, e.g., to spiral-wound membranes, or even for different applications, such as soil mechanics and petroleum engineering [52].

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript:

CFD
Computational fluid dynamics SHFM Submerged hollow fiber membrane