Abstract
This paper presents an analytical model of substrate mass transfer through the lumen of a membrane bioreactor. The model is a solution of the convective-diffusion equation in two dimensions using a regular perturbation technique. The analysis accounts for radial-convective flow as well as axial diffusion of the substrate specie. The model is applicable to the different modes of operation of membrane bioreactor (MBR) systems (e.g., dead-end, open-shell, or closed-shell mode), as well as the vertical or horizontal orientation. The first-order limit of the Michaelis-Menten equation for substrate consumption was used to test the developed model against available analytical results. The results obtained from the application of this model, along with a biofilm growth kinetic model, will be useful in the derivation of an efficiency expression for enzyme production in an MBR.
1. Introduction
Since the first uses of hollow-fiber membrane bioreactors (MBRs) to immobilize whole cells were reported in the early 1970s, this technology has been used in as wide ranging applications as enzyme production to bone tissue engineering. One of the current research areas of interest into biofilm-attached membrane bioreactors (MBRs) is the development of cost-effective and environmentally friendly methods of producing various primary and secondary metabolites from bacterial, fungal, and yeast cells. These include: manganese and lignin peroxidase, secreted by the fungus Phanerochaete chrysosporium [1, 2]; actinorhodin, a noncommercial antibiotic produced by the filamentous bacterium Streptomyces coelicolor [3]; glutamic acid, an ingredient in flavour enhancers of meats and vegetables, secreted by the bacterium Corynebacterium glutamicum [4]; ethanol, extracted from the yeast Saccharomyces cerevisiae [5]; and many others. With the exception of ethanol, these bioproducts are generally classified as products of intermediate value [6]. It has been reported that bioreactor productivity, in the production of these types of products, greatly impacts on the product cost [7].
The productivity of biofilm-attached MBRs is determined in large by the biomass growth, and one of the most important factors that influence biomass growth is the availability and transport of nutrients through the bioreactor [8, 9]. The momentum transfer of solutes through MBRs has been thoroughly studied, from a theoretical and experimental perspective, for a number of configurations [10–15]. Similarly, the mass transfer has received considerable attention [8, 9, 16–22]. With the exception of the models developed by Heath and Belfort [17]; Li and Tan [19]; and Willaert et al. [22], the mass transfer models were solved using numerical procedures such as finite difference schemes and control volumes. A difficulty in implementing such schemes is the choice of the appropriate technique for a specific MBR system [21], and these techniques are subject to discretization errors and stringent stability criterion. In the models presented by Heath and Belfort [17]; Willaert et al. [22]; and Li and Tan [19], the convective-diffusion equation governing mass transfer was solved analytically. These authors, however, neglected the effects of axial diffusion and radial convective flow in their models. Both these assumptions may not be justified in all cases. A number of theoretical and experimental investigations have demonstrated the significance of radial convective flows in improving MBR efficiencies [9, 15, 21]. In the dead-end ultrafiltration mode, particularly, the assumption of negligible radial convective flow is not justifiable. At axial Peclet numbers (Peu) smaller than unity large concentration gradients exist, and under these circumstances ignoring axial diffusion is also not justified.
The current study presents an analytical solution of the convective-diffusion equation, for solute transport trough a single fiber isotropic capillary membrane, in two dimensions. This study will not include the growth kinetics of the microorganism, as conversion is assumed to take place in the shell side of the MBR; the current analysis is restricted to the lumen side. For comparison with literature models, however, the first-order limit of the Michaelis-Menten equation will be superimposed on the developed model in the results section.
2. Model Development
2.1. The Membrane Gradostat Reactor
The models developed in this study are applicable to a hollow-fiber MBR system, consisting of either a single fibre or a bundle of fibres; with nutrient flowing on the lumen side of the membrane and the micro-organism immobilised either on the lumen side or on the shell side. The notation used, however, is specifically for a single hollow fiber membrane gradostat reactor (MGR). The construction of the MGR, as patented by Edwards et al. [23], is illustrated schematically in Figure 1. It consists of a single hollow-fibre, made of surface modified polysulphone, encased in a glass bioreactor. The membranes are asymmetric and characterized by an internally skinned and externally unskinned region of microvoids; approximately 0.15 mm long and 0.015 mm thick. These membranes have inner and outer diameters of approximately 1.395 mm and 1.925 mm, respectively. The nutrient solution permeates from the lumen side to the shell side of the MGR due to the transmembrane pressure gradient. The micro-organism is immobilised on the shell side of the MGR. Humidified air is supplied on the shell side, and two pressure transducers are fitted at the inlet and outlet of the MGR as shown in Figure 1.
2.2. Model Assumptions
The theoretical models to be developed will be based on the following conditions of operation and assumptions: () the system is isothermal, meaning the energy equation has been decoupled from the mass and momentum transfer; () the flow regime within the membrane lumen is fully developed, laminar, homogeneous, and at steady state; () the physical and transport parameters (e.g., density, viscosity, and diffusivity) are constant; () in the dense and spongy layers of the membrane matrix the flow is only one dimensional (i.e., there are no axial components of the velocity profiles in the membrane matrix); and () the aspect ratio of the membrane is much smaller than unity. The aspect ratio, , is the ratio of the membrane inner radius to the effective membrane length (i.e., ), and if it is much smaller than unity then normal stress effects are negligible in the momentum transfer analysis.
3. Mathematical Formulation
The starting point of the analysis is the convective-diffusion equation [24]: where is the local substrate concentration; is time; is the substrate diffusivity, assumed to be constant; and , the rate of substrate production (or consumption), is a function of the local biofilm density, and time. Equation (1), for steady state, two-dimensional flow, without reaction, in cylindrical co ordinates may be written as It is convenient to express this equation in dimensionless form by introducing the following dimensionless variables: The expressions of and in (3) are obtained from the momentum transfer analysis given in Appendix A. Substituting the dimensionless variables in (3) into (2) results in where the axial and radial Peclet numbers () are defined as The boundary conditions which match the imposed operating conditions of the MBR system are presented in Table 1. Boundary condition 1 (B.C.1) corresponds to a uniform inlet substrate concentration; B.C.2 and B.C.8 corresponds to cylindrical symmetry at the centre of the membrane lumen; B.C.3 indicates that the substrate concentration is only a function of the axial spatial coordinate at the membrane wall; B.C.4 corresponds to continuity of the substrate flux at the lumen-matrix interface; B.C.5 indicates a distance along the length of the membrane above which the axial concentration gradient becomes zero. The boundary conditions B.C.7–11 are employed in the solution of the velocity profiles given in Appendix A. A cross section of the MGR illustrating the three regions of the reactor, with the notation used in the model development, is shown in Figure 2.
If in (4) is radially averaged, to become , then the left-hand side (LHS) of (4) is only a function of Z and the right-hand side (RHS) only a function of . This can only be true if both the LHS and RHS are independent of the variables and . Equation (4) may therefore be solved by separation of variables to give a solution of the form Substituting (6) into (4) gives The equating of the two ordinary differential equations (ODEs) to the arbitrary constant in (7) is due to the fact that the two ODEs are independent of the variables and .
3.1. Solution of the Axial Concentration Function
To solve for the axial function of the substrate concentration profile, the ODE on the LHS of (7) is considered The radially-averaged axial velocity in (8) is defined as where and are the Froude and Reynolds number, respectively, given by where is the gravitational acceleration; is the solution density; and is the solution dynamic viscosity. The solution of the axial velocity in (9) is given in Appendix A as with the axial pressure gradient given by where is the dimensionless hydrostatic pressure; is the dimensionless transmembrane pressure; is the dimensionless entrance pressure drop; and is the dimensionless membrane hydraulic permeability. The entrance pressure drop in (12) is given by where is the fraction retentate ( for the dead-end mode and for the closed-shell mode). The membrane hydraulic permeability in (12) is much smaller than unity (), therefore, this equation can be approximated by the following expression: This approximation makes (8) a confluent hypergeometric type differential equation. This is more evident if the following sequential substitutions are made.
Substitution 1
The substitution in (15) represents the axial gradient (or driving force) of the substrate concentration profile. The axial distance along the membrane length at which this gradient is zero represents the point at which the axial function of the concentration profile, , is constant. This axial distance, , is obtained by equating (15) to zero as follows:
The substitution in (15) transforms (8) to
where
Substitution 2
This substitution transforms (17) to
Equation (20) is the standard Kummer hypergeometric equation and has two solutions, the Kummer function of the first kind and the Tricomi function , respectively [25]:
where
where is the gamma function. Therefore, the solution of (20) becomes
where and are coefficients obtained from the inlet condition ((6a) in Table 1). The Tricomi function approaches infinity as values of approach zero [26], therefore, the coefficient in (23) must be zero for this equation to satisfy ((6a) in Table 1). The coefficient is obtained from the inlet condition ((6a) in Table 1):
where
From the definition in (21), the Kummer function is equal to 1 for all real values of and (where ). The two piecewise solutions of the axial function of the dimensionless concentration profile, (), are therefore
3.2. Solution of the Radial Concentration Function
3.2.1. Zero-Order Approximation
To solve for the radial function of the substrate concentration profile, the ODE on the RHS of (7) is solved The radial velocity in (27) is given in Appendix A as From the approximation in(14) Equation (27) is solved by a regular perturbation technique: The magnitude of the membrane hydraulic permeability is very small; hence, validity of the perturbation method is assured. The equations to solve for the zero-order approximation, , of (27) are Equation (31) is Bessel’s differential equation and has a standard solution of the form As approaches zero in (33), the function tends to minus infinity; and therefore, must be zero for the equation to satisfy B.C.6 at is equal to zero. Imposing B.C.3, the dimensionless concentration is only a function of the axial coordinate at ; therefore, the function should be equal to the constant at , but is an oscillating function that can have a number of roots that satisfy the condition of B.C.3: The solution for the coefficient in (35) is given in Appendix B.1; the eigenvalues are derived from B.C.4 and are given in Appendix B.2.
3.3. First-Order and Second-Order Approximations
The equations to solve for the first-order approximation, , of (27) are where Equation (36) is evaluated by making use of the following identity of Bessel functions [25]: Substituting (39) into (36) results in the following inhomogeneous O.D.E: Equation (40) is further simplified by making use of the following substitution: This substitution simplifies (40) to Some mathematical architecture is required to solve (42) and this is described in Appendix C.1. The solution of this equation is given in Appendix C.2 as where The differential equations to solve for the second-order approximation, , of (27) are The mathematical architecture required for the solution of these equations is also given in Appendix C.1, and the solution is given in Appendix C.3: Substituting (26) and (30) into (6), the solution of the dimensionless luminal concentration profile is therefore The perturbation solution is only extended up to second-order approximation; therefore, (48) reduces to
4. Results
4.1. Model Validation
The small diameters of the capillary membranes used in the construction of the MBRs render it a very difficult task to validate the accuracy of the developed models experimentally, and thus the predictive power of the models will be compared to currently available results for similar MBR systems; in particular the model developed by Heath and Belfort [17]. The approach used by these authors in solving the convective-diffusion equation was to assume a constant cross-sectional concentration profile in the membrane lumen and matrix regions, and to solve the simplified form of the equation. The lumen side axial concentration profile was then obtained from a mass flux balance. This approach is equivalent to assuming that the substrate consumption takes place in the lumen side of the MBR; in such case the convective-diffusion equation given in (1) becomes where is the maximum rate of reaction and is the Michaelis constant. This equation is made dimensionless by introducing the variables in (3) and (5), and two-additional variables the Thiele modulus and the dimensionless Michaelis constant , respectively: Equation (50) then becomes Assuming the first-order limit of the Michaelis-Menten equation (i.e., ) the solution of (52) takes the same general approach as that of (4). The solution of this equation is identical to (49) with an adjustment of the axial function and the coefficient of the perturbation solution to account for the reaction rate: where The input variables used in validating (53) are obtained from Heath and Belfort [17] and are given in Table 2. The values of the membrane hydraulic permeability (), the fraction retentate (), the lumen-side inlet hydrostatic pressure (), the shell-side hydrostatic pressure (), the solution density (), and viscosity () are not specified by these authors; therefore, typical operational values of these parameters will be used. A transmembrane (TMP) pressure of 5 kPa is assumed across the MBR, and the solution properties are assumed to be those of water at 30°C.
In Figure 3, it can be seen that (53) compares satisfactorily with the model of Heath and Belfort [17] for the parameter values listed in Table 2. At the centre of the MBR () Equation (53) predicts that the concentration decreases axially to 34% of its original value when the applied TMP is 5kPa, and at the membrane wall () the concentration decreases up to 24% of its original value as shown in Figure 4. The model of Heath and Belfort [17] predicts a 44% decrease for all values of for the corresponding conditions. It is important to note, however, that only the resulting trends from the two models can be compared, since (53) requires a more detailed description of the MBR system than that required for the model of Heath and Belfort [17].
5. Conclusion
A rigour analytical mathematical model for substrate concentration profiles in the lumen of a hollow fiber MBR was developed. The model was based on the solution of the convective-diffusion equation in dimensionless form. The model allows evaluation of the influence of the general operating parameters of a MBR on the concentration profiles. These parameters are the fraction retentate (); the membrane hydraulic permeability (); the axial and radial Peclet numbers (); the Thiele modulus (); the fluid properties; and the dimensions of the MBR. The developed model can be further used to evaluate reactor performance from basic principles, since it allows analytical evaluation of performance parameters (e.g., the performance index and the effectiveness factor). The current paper is a precursor to a paper by the same authors on bioreactor performance. The models developed in the current paper will be used to develop expressions for the effectiveness factor and performance index (PI) of a capillary membrane gradostat reactor.
Appendices
A. Momentum Transfer Analysis
A.1. Momentum Transfer Analysis
The z-component of the nonconservative form of the Navier-Stokes equation in cylindrical coordinates is made dimensionless by introducing the variables in (4) and the following additional variables: where and are the hydrostatic and dimensionless hydrostatic pressures, respectively; and are the time and dimensionless time, respectively; is the solution dynamic viscosity; and is the Froude number. The dimensionless form of the Navier-Stokes thus becomes and the continuity equation Equation (A.2) is solved following the procedure of Godongwana et al. [13], and the boundary conditions listed in Table 1 (B.C.7 and B.C.8). Ignoring normal stresses and considering steady-state conditions, the solution of (A.2) is given by The dimensionless radial velocity profile is obtained by substituting (A.4) into (A.3) and imposing B.C.9: The dimensionless pressure profile is obtained by imposing B.C.10, where the matrix velocity is governed by Darcy’s law: where the dimensionless membrane hydraulic permeability is given by
Substituting (A.5) and (A.6) into (6j) results in Equation (A.8) is solved by applying B.C.11 to give: where the dimensionless entrance pressure drop in (A.9) is obtained by defining a fraction retentate, , as the ratio of the exit to the inlet velocity, that is, :
B. Solution of Constants
B.1. Solution of the Coefficient
Imposing the inlet condition (6a) into (48) gives For the zero-order approximation of the radial function , (B.1) becomes To solve for in (B.2) both the LHS and the RHS are multiplied by and integrated with respect to over the interval 0-1, The RHS of (B.3) is evaluated by making use of the following property of Bessel functions: The R.H.S of (B.3) then becomes The L.H.S of (B.3) may be evaluated by making use of Lommel integrals: To give the following equation: Substituting (B.5) and (B.7) back into the RHS and LHS of (B.3), respectively
B.2. Solution of the Eigenvalues
The solution of the eigenvalues is obtained from B.C.4 in Table 1. When the enzyme is immobilized on the lumen side of the membrane this equation reduces to Substituting the expression for C, (49) into (B.9) and limiting the perturbation to first-order approximation: The eigenvalues of (B.8) are values of that satisfy (B.10). The first 10 of these values are listed in Table 3 for the parameter values in Table 2.
C. Laplace Transform Solutions
C.1. Mathematical Architecture for the Solution of
We define a parameter in Laplace space as The following properties of are all noteworthy, differentiation: Integration: Inverse Laplace Transform: Laplace Transforms involving the First-Order Bessel function: More relations involving : We also define the following polynomial in , which is significant for the first- and second-order approximation of in Sections C.2 and C.3, respectively, Then it is easy to show that the second derivative of the product with respect to the Laplace variable is The fourth derivative of the product with respect to the Laplace variable is Both (C.9) and (C.10) will become important in Appendix C.3.
C.2. Solution of the First-Order Approximation Function, , in (42)
If the function is taken as the Laplace transform of the function , that is, , then the Laplace transform of (42) yields In terms of the parameter , defined in Appendix C.1, this equation may be written as The right-hand side of (C.12) is evaluated from the relations involving in Appendix C.1 (C.7). Multiplying through by (C.12) becomes Integrating (C.13) with respect to , to find the Laplace space solution: For convenience of expressing the solutions of the first- and second-order approximations, and , respectively, the following constants and function are defined:where The polynomial in (C.15c) was defined in (C.8). The inverse Laplace transform of the function in (C.14) is the solution of :
C.3. Solution of the Second-Order Approximation Function, , in (45)
Similar to the first-order approximation, the Laplace transform of the second-order approximation, for , yields In terms of the parameter , defined in Appendix C.1, this equation may be written as The expressions for the fourth derivative and second derivative of are given in Section C.1, (C.9)-(C.10). The first-order differential equation in has an integrating factor , so again (C.19) is multiplied by to obtain The integral of (C.20), after simplifying by grouping terms with like powers of , is Recalling that , the second-order approximation of (27) is simply the inverse Laplace transform of (C.21):
Nomenclature
| : | Dimensionless entrance pressure drop |
| : | Constants of integration of Bessel’s equation, , 2 |
| : | Substrate concentration (g dm−3) |
| : | Substrate feed concentration (g dm−3) |
| : | Dimensionless substrate concentration |
| : | Substrate diffusivity (m2 s−1) |
| : | Fraction retentate |
| : | Coefficients of the solution of Kummer’s confluent hypergeometric equation, , 2 |
| : | Froude number |
| : | Dimensionless axial concentration function |
| : | Gravitational acceleration (m s−2) |
| : | Laplace transform of the first-order approximation of the function |
| : | Laplace transform of the second-order approximation of the function |
| : | Constants in the first and second-order approximations of the function , , 2, 3 |
| : | Bessel function of order of the first kind |
| : | Mass transfer coefficient (m s−1) |
| : | Membrane hydraulic permeability (m Pa−1 s−1) |
| : | Michaelis constant (g dm−3) |
| : | Dimensionless Michaelis constant |
| : | Membrane effective length (m) |
| : | Kummer function of the first kind |
| : | Hydrostatic pressure (Pa) |
| : | Dimensionless hydrostatic pressure |
| : | Axial Peclet number |
| : | Radial Peclet number |
| : | Rate of substrate production/consumption (g dm−3 s−1) |
| : | Radial spatial coordinate (m) |
| : | Dimensionless radial spatial coordinate |
| : | Membrane lumen radius (m) |
| : | Reynolds number |
| : | Sherwood number |
| : | Time (s) |
| : | Dimensionless radial concentration function |
| : | Axial velocity (m s−1) |
| : | Feed axial velocity (m s−1) |
| : | Dimensionless axial velocity |
| : | Radial velocity (m s−1) |
| : | Permeation velocity (m s−1) |
| : | Dimensionless radial velocity |
| : | Maximum rate of reaction (g dm−3 s−1) |
| : | Substitution variable |
| : | Bessel function of order of the second kind |
| : | Axial spatial coordinate (m) |
| : | Dimensionless axial spatial coordinate |
| : | Dimensionless axial distance at which the concentration gradient is zero. |
| : | First parameter in the Kummer functions of the first and second kind |
| : | Dimensionless transmembrane pressure |
| : | Lumped parameter in(38) |
| : | Membrane porosity |
| : | Thiele modulus |
| : | Tricomi function/Kummer function of the second kind |
| : | Second parameter in the Kummer functions of the first and second kind |
| : | Gamma function, |
| : | Aspect ratio |
| : | Dimensionless membrane hydraulic permeability |
| : | Eigen values, |
| : | Solution dynamic viscosity (Pa s) |
| : | Substitution variable |
| : | Solution density (kg m−3) |
| : | Dimensionless time |
| : | Axial gradient/driving force of substrate concentration profile. |
| 0: | Membrane entrance |
| 1: | Membrane exit |
| : | Shell-side of the MBR. |
Acknowledgments
The authors would like to thank the Deutscher Akademischer Austausch Dienst (DAAD) for financial support, Dr. Aletta van der Merwe for her contribution in validating some of the mathematical principles.