# A Solution of the Convective-Diffusion Equation for Solute Mass Transfer inside a Capillary Membrane Bioreactor

Authors: B. Godongwana; D. Solomons; M. S. Sheldon
Journal: International Journal of Chemical Engineering
 (2010)

---

## 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.


---

## Body

## 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 fungusPhanerochaete 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.Figure 1
A schematic diagram of the single capillary membrane gradostat reactor (MGR).


### 2.2. Model Assumptions

The theoretical models to be developed will be based on the following conditions of operation and assumptions: (1) the system is isothermal, meaning the energy equation has been decoupled from the mass and momentum transfer; (2) the flow regime within the membrane lumen is fully developed, laminar, homogeneous, and at steady state; (3) the physical and transport parameters (e.g., density, viscosity, and diffusivity) are constant; (4) 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 (5) 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., RL/L), and if it is much smaller than unity then normal stress effects are negligible in the momentum transfer analysis.


## 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.Figure 1
A schematic diagram of the single capillary membrane gradostat reactor (MGR).


## 2.2. Model Assumptions

The theoretical models to be developed will be based on the following conditions of operation and assumptions: (1) the system is isothermal, meaning the energy equation has been decoupled from the mass and momentum transfer; (2) the flow regime within the membrane lumen is fully developed, laminar, homogeneous, and at steady state; (3) the physical and transport parameters (e.g., density, viscosity, and diffusivity) are constant; (4) 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 (5) 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., RL/L), 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]: (1)DcDt=DAB∇2c+rA,
where c is the local substrate concentration; t is time; DAB is the substrate diffusivity, assumed to be constant; and rA, 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(2)u∂c∂z-DAB∂2c∂z2=DABr(∂c∂r+r∂2c∂r2)-v∂c∂r.
It is convenient to express this equation in dimensionless form by introducing the following dimensionless variables:(3)U=uu0,V=vv0,C=cc0,Z=zL,R=rRL,φ=RLL.
The expressions of U and V in (3) are obtained from the momentum transfer analysis given in Appendix A. Substituting the dimensionless variables in (3) into (2) results in(4)φPeuU∂C∂Z-φ2∂2C∂Z2=1R(∂C∂R+R∂2C∂R2)-PevV∂C∂R,
where the axial and radial Peclet numbers (Peu,v) are defined as(5)Peu=u0RLDAB,Pev=v0RLDAB.
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.Table 1
The boundary conditions of the MBR.
B.C.R, ZC, U, V, PLRangeEquationB.C.1Z=0C=10≤R≤ 1(6a)B.C. 2R=0∂C∂R=00≤Z≤ 1(6b)B.C. 3R=1C=C(Z)0≤Z≤ 1(6c)B.C. 4R=1-∂C∂R+PevVC=Sh[C(Z,1)-C(Z,R2)]0≤Z≤ 1(6d)B.C. 5Z≥φ16βκ(ReFr-aφ)C(Z) = constant0≤R≤ 1(6e)B.C. 6R=0C = finite0≤Z≤ 1(6f)B.C. 7R=1U=00≤Z≤ 1(6g)B.C. 8R=0∂U∂R=00≤Z≤ 1(6h)B.C. 9R=0V=00≤Z≤ 1(6i)B.C.10R=1V=VM0≤Z≤ 1(6j)B.C.11Z=0P=P0; P′=φa0≤R≤ 1(6k)Figure 2
A cross section of the single fibre MGR.IfU in (4) is radially averaged, to become Uav, then the left-hand side (LHS) of (4) is only a function of Z and the right-hand side (RHS) only a function of R. This can only be true if both the LHS and RHS are independent of the variables R and Z. Equation (4) may therefore be solved by separation of variables to give a solution of the form(6)C=F(Z)T(R).
Substituting (6) into (4) gives(7)φPeuUavFdFdZ-φ2Fd2FdZ2=1RT(dTdR+Rd2TdR2)-PevVTdTdR=-λ2.
The equating of the two ordinary differential equations (ODEs) to the arbitrary constant -λ2 in (7) is due to the fact that the two ODEs are independent of the variables R and Z.


### 3.1. Solution of the Axial Concentration FunctionF(Z)

To solve for the axial functionF(Z) of the substrate concentration profile, the ODE on the LHS of (7) is considered(8)d2FdZ2-PeuUavφdFdZ-λ2φ2F=0.
The radially-averaged axial velocity Uav in (8) is defined as(9)Uav=2∫01URdR=-18(dPdZ-ReFr),
where Fr and Re are the Froude and Reynolds number, respectively, given by(10)Fr=u02gRL,Re=ρu0RLμ,
where g is the gravitational acceleration; ρ is the solution density; and μ is the solution dynamic viscosity. The solution of the axial velocity U in (9) is given in Appendix A as(11)U=-14(1-R2)(dPdZ-ReFr),
with the axial pressure gradient given by(12)dPdZ=4φ-1κβsinh(4φ-1κ)Z+φacosh(4φ-1κ)Z,
where P is the dimensionless hydrostatic pressure; β is the dimensionless transmembrane pressure; a is the dimensionless entrance pressure drop; and κ is the dimensionless membrane hydraulic permeability. The entrance pressure drop a in (12) is given by(13)a=4φ-1κβsinh(4φ-1κ)-ReFr-1(1-f)φ[f-cosh(4φ-1κ)],
where f is the fraction retentate (f=0 for the dead-end mode and f=1 for the closed-shell mode). The membrane hydraulic permeability κ in (12) is much smaller than unity (κ≪1), therefore, this equation can be approximated by the following expression:(14)dPdZ≈φa+16βφ-1κZ.
This approximation makes (8) a confluent hypergeometric type differential equation. This is more evident if the following sequential substitutions are made.Substitution 1
(15)ξ={-Peu8φ(aφ+16βκZφ-ReFr),0≤Z<Z0,0,Z0≤Z≤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, F(Z), is constant. This axial distance, Z0, is obtained by equating (15) to zero as follows:
(16)∴ξ=0atZ0=φ16βκ(ReFr-φa).
The substitution in (15) transforms (8) to
(17)d2Fdξ2-AξdFdξ-A2λ2φ2F=0,
where
(18)A=-φ2(2Peuβκ)-1.Substitution 2
(19)θ=12Aξ2.
This substitution transforms (17) to
(20)θd2Fdθ2+(12-θ)dFdθ-Aλ22φ2F=0.
Equation (20) is the standard Kummer hypergeometric equation and has two solutions, the Kummer function of the first kind M(α,γ,θ) and the Tricomi function Φ(α,γ,θ), respectively [25]:
(21)M(α,γ,θ)=1+αγθ+(α)2(γ)22!θ2+⋯+(α)n(γ)nn!θn,
where
(22)(α)n=α(α+1)(α+2)+⋯(α+n-1),(α)0=1,Φ(α,γ,θ)=πsinπγ{M(α,γ,θ)Γ(1+α-γ)Γ(γ)-θ1-γM(1+α-γ,2-γ,θ)Γ(α)Γ(2-γ)},
where Γ(n) is the gamma function. Therefore, the solution of (20) becomes
(23)F(θ)=F0M(Aλ22φ2,12,θ)+F1Φ(Aλ22φ2,12,θ),
where F0and F1 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 F1 in (23) must be zero for this equation to satisfy ((6a) in Table 1). The coefficient F0 is obtained from the inlet condition ((6a) in Table 1):
(24)C=F0M(Aλ22φ2,12,θ0)T(R)=1atZ=0,
where
(25)θ0=-Peu256βκ(aφ-ReFr)2.
From the definition in (21), the Kummer function M(α,γ,0) is equal to 1 for all real values of α and γ (where γ≠0). The two piecewise solutions of the axial function of the dimensionless concentration profile, F(θ), are therefore
(26)F(θ)={F0M(Aλ22φ2,12,θ),0≤Z<Z0,F0,Z0≤Z≤1.


### 3.2. Solution of the Radial Concentration FunctionT(R)

#### 3.2.1. Zero-Order Approximation

To solve for the radial functionT(R) of the substrate concentration profile, the ODE on the RHS of (7) is solved(27)d2TdR2+(1R-PevV)dTdR+λ2T=0.
The radial velocity V in (27) is given in Appendix A as(28)V=φ(u0v0)[R8(1-R22)]d2PdZ2.
From the approximation in(14)(29)d2PdZ2≈16βφ-1κ,κ≪1.
Equation (27) is solved by a regular perturbation technique:(30)T(R)=∑n=0∞κnTn.
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, T0, of (27) are(31)d2T0dR2+1RdT0dR+λ2T0=0;(32)T0(0)-B1=0;dT0(0)dR-B2=0
Equation (31) is Bessel’s differential equation and has a standard solution of the form(33)T0=B1J0(λR)+B2Y0(λR).
As R approaches zero in (33), the function Y0 tends to minus infinity; and therefore, B2 must be zero for the equation to satisfy B.C.6 at R is equal to zero.(34)∴T0=B1J0(λR).
Imposing B.C.3, the dimensionless concentration C is only a function of the axial coordinate Z at R=1; therefore, the function T0 should be equal to the constant B1J0(λ) at R=1, but J0(λ) is an oscillating function that can have a number of roots that satisfy the condition of B.C.3:(35)∴T0=∑m=1∞B1mJ0(λmR).
The solution for the coefficient B1m in (35) is given in Appendix B.1; the eigenvalues λm 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,T1, of (27) are(36)d2T1dR2+1RdT1dR+λm2T1=δ[R(1-R22)]dJ0(λmR)dR,(37)T1(0)=0,dT1(0)dR=0,
where(38)δ=2Pevβ(u0v0)∑m=1∞B1m.
Equation (36) is evaluated by making use of the following identity of Bessel functions [25]:(39)dJ0(λmR)dR=-λmJ1(λmR).
Substituting (39) into (36) results in the following inhomogeneous O.D.E:(40)d2T1dR2+1RdT1dR+λm2T1=-δλm[R(1-R22)]J1(λmR).
Equation (40) is further simplified by making use of the following substitution:(41)x=λmR.
This substitution simplifies (40) to(42)xd2T1dx2+dT1dx+xT1=-δx2λm2(1-x22λm2)J1(x).
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(43)T1(x)=i3[x2J2(x)3!!+i1x3J3(x)5!!+i2x4J4(x)7!!],
where(44)i1=-203λm2,i2=354λm2,i3=-3δ4λm2.
The differential equations to solve for the second-order approximation, T2, of (27) are(45)xd2T2dx2+dT2dx+xT2=-δx2λm2(1-x22λm2)T1′(x),(46)T2(0)=0,dT2(0)dx=0.
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:(47)T2(x)=(3δ22λm6){203λm2x3J3(x)5!!-[105-73λm2(6i1-5)]x4J4(x)7!!+910[7(60-48i1)+2λm2(8i2-7i1)]x5J5(x)9!!-9912[(35-112i1+80i2)+2λm2i2]x6J6(x)11!!-9·11·1314(7i1-20i2)x7J7(x)13!!-9·11·13·1516·i2x8J8(x)15!!}.
Substituting (26) and (30) into (6), the solution of the dimensionless luminal concentration profile is therefore(48)C(θ,x)=∑m=1∞∑n=0∞Fm(θ)Tmn(x)κn.
The perturbation solution is only extended up to second-order approximation; therefore, (48) reduces to(49)C(θ,x)=∑m=1∞Fm(θ)[Tm0(x)κ0+Tm1(x)κ1+Tm2(x)κ2].


## 3.1. Solution of the Axial Concentration FunctionF(Z)

To solve for the axial functionF(Z) of the substrate concentration profile, the ODE on the LHS of (7) is considered(8)d2FdZ2-PeuUavφdFdZ-λ2φ2F=0.
The radially-averaged axial velocity Uav in (8) is defined as(9)Uav=2∫01URdR=-18(dPdZ-ReFr),
where Fr and Re are the Froude and Reynolds number, respectively, given by(10)Fr=u02gRL,Re=ρu0RLμ,
where g is the gravitational acceleration; ρ is the solution density; and μ is the solution dynamic viscosity. The solution of the axial velocity U in (9) is given in Appendix A as(11)U=-14(1-R2)(dPdZ-ReFr),
with the axial pressure gradient given by(12)dPdZ=4φ-1κβsinh(4φ-1κ)Z+φacosh(4φ-1κ)Z,
where P is the dimensionless hydrostatic pressure; β is the dimensionless transmembrane pressure; a is the dimensionless entrance pressure drop; and κ is the dimensionless membrane hydraulic permeability. The entrance pressure drop a in (12) is given by(13)a=4φ-1κβsinh(4φ-1κ)-ReFr-1(1-f)φ[f-cosh(4φ-1κ)],
where f is the fraction retentate (f=0 for the dead-end mode and f=1 for the closed-shell mode). The membrane hydraulic permeability κ in (12) is much smaller than unity (κ≪1), therefore, this equation can be approximated by the following expression:(14)dPdZ≈φa+16βφ-1κZ.
This approximation makes (8) a confluent hypergeometric type differential equation. This is more evident if the following sequential substitutions are made.Substitution 1
(15)ξ={-Peu8φ(aφ+16βκZφ-ReFr),0≤Z<Z0,0,Z0≤Z≤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, F(Z), is constant. This axial distance, Z0, is obtained by equating (15) to zero as follows:
(16)∴ξ=0atZ0=φ16βκ(ReFr-φa).
The substitution in (15) transforms (8) to
(17)d2Fdξ2-AξdFdξ-A2λ2φ2F=0,
where
(18)A=-φ2(2Peuβκ)-1.Substitution 2
(19)θ=12Aξ2.
This substitution transforms (17) to
(20)θd2Fdθ2+(12-θ)dFdθ-Aλ22φ2F=0.
Equation (20) is the standard Kummer hypergeometric equation and has two solutions, the Kummer function of the first kind M(α,γ,θ) and the Tricomi function Φ(α,γ,θ), respectively [25]:
(21)M(α,γ,θ)=1+αγθ+(α)2(γ)22!θ2+⋯+(α)n(γ)nn!θn,
where
(22)(α)n=α(α+1)(α+2)+⋯(α+n-1),(α)0=1,Φ(α,γ,θ)=πsinπγ{M(α,γ,θ)Γ(1+α-γ)Γ(γ)-θ1-γM(1+α-γ,2-γ,θ)Γ(α)Γ(2-γ)},
where Γ(n) is the gamma function. Therefore, the solution of (20) becomes
(23)F(θ)=F0M(Aλ22φ2,12,θ)+F1Φ(Aλ22φ2,12,θ),
where F0and F1 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 F1 in (23) must be zero for this equation to satisfy ((6a) in Table 1). The coefficient F0 is obtained from the inlet condition ((6a) in Table 1):
(24)C=F0M(Aλ22φ2,12,θ0)T(R)=1atZ=0,
where
(25)θ0=-Peu256βκ(aφ-ReFr)2.
From the definition in (21), the Kummer function M(α,γ,0) is equal to 1 for all real values of α and γ (where γ≠0). The two piecewise solutions of the axial function of the dimensionless concentration profile, F(θ), are therefore
(26)F(θ)={F0M(Aλ22φ2,12,θ),0≤Z<Z0,F0,Z0≤Z≤1.


## 3.2. Solution of the Radial Concentration FunctionT(R)

### 3.2.1. Zero-Order Approximation

To solve for the radial functionT(R) of the substrate concentration profile, the ODE on the RHS of (7) is solved(27)d2TdR2+(1R-PevV)dTdR+λ2T=0.
The radial velocity V in (27) is given in Appendix A as(28)V=φ(u0v0)[R8(1-R22)]d2PdZ2.
From the approximation in(14)(29)d2PdZ2≈16βφ-1κ,κ≪1.
Equation (27) is solved by a regular perturbation technique:(30)T(R)=∑n=0∞κnTn.
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, T0, of (27) are(31)d2T0dR2+1RdT0dR+λ2T0=0;(32)T0(0)-B1=0;dT0(0)dR-B2=0
Equation (31) is Bessel’s differential equation and has a standard solution of the form(33)T0=B1J0(λR)+B2Y0(λR).
As R approaches zero in (33), the function Y0 tends to minus infinity; and therefore, B2 must be zero for the equation to satisfy B.C.6 at R is equal to zero.(34)∴T0=B1J0(λR).
Imposing B.C.3, the dimensionless concentration C is only a function of the axial coordinate Z at R=1; therefore, the function T0 should be equal to the constant B1J0(λ) at R=1, but J0(λ) is an oscillating function that can have a number of roots that satisfy the condition of B.C.3:(35)∴T0=∑m=1∞B1mJ0(λmR).
The solution for the coefficient B1m in (35) is given in Appendix B.1; the eigenvalues λm are derived from B.C.4 and are given in Appendix B.2.


## 3.2.1. Zero-Order Approximation

To solve for the radial functionT(R) of the substrate concentration profile, the ODE on the RHS of (7) is solved(27)d2TdR2+(1R-PevV)dTdR+λ2T=0.
The radial velocity V in (27) is given in Appendix A as(28)V=φ(u0v0)[R8(1-R22)]d2PdZ2.
From the approximation in(14)(29)d2PdZ2≈16βφ-1κ,κ≪1.
Equation (27) is solved by a regular perturbation technique:(30)T(R)=∑n=0∞κnTn.
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, T0, of (27) are(31)d2T0dR2+1RdT0dR+λ2T0=0;(32)T0(0)-B1=0;dT0(0)dR-B2=0
Equation (31) is Bessel’s differential equation and has a standard solution of the form(33)T0=B1J0(λR)+B2Y0(λR).
As R approaches zero in (33), the function Y0 tends to minus infinity; and therefore, B2 must be zero for the equation to satisfy B.C.6 at R is equal to zero.(34)∴T0=B1J0(λR).
Imposing B.C.3, the dimensionless concentration C is only a function of the axial coordinate Z at R=1; therefore, the function T0 should be equal to the constant B1J0(λ) at R=1, but J0(λ) is an oscillating function that can have a number of roots that satisfy the condition of B.C.3:(35)∴T0=∑m=1∞B1mJ0(λmR).
The solution for the coefficient B1m in (35) is given in Appendix B.1; the eigenvalues λm 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,T1, of (27) are(36)d2T1dR2+1RdT1dR+λm2T1=δ[R(1-R22)]dJ0(λmR)dR,(37)T1(0)=0,dT1(0)dR=0,
where(38)δ=2Pevβ(u0v0)∑m=1∞B1m.
Equation (36) is evaluated by making use of the following identity of Bessel functions [25]:(39)dJ0(λmR)dR=-λmJ1(λmR).
Substituting (39) into (36) results in the following inhomogeneous O.D.E:(40)d2T1dR2+1RdT1dR+λm2T1=-δλm[R(1-R22)]J1(λmR).
Equation (40) is further simplified by making use of the following substitution:(41)x=λmR.
This substitution simplifies (40) to(42)xd2T1dx2+dT1dx+xT1=-δx2λm2(1-x22λm2)J1(x).
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(43)T1(x)=i3[x2J2(x)3!!+i1x3J3(x)5!!+i2x4J4(x)7!!],
where(44)i1=-203λm2,i2=354λm2,i3=-3δ4λm2.
The differential equations to solve for the second-order approximation, T2, of (27) are(45)xd2T2dx2+dT2dx+xT2=-δx2λm2(1-x22λm2)T1′(x),(46)T2(0)=0,dT2(0)dx=0.
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:(47)T2(x)=(3δ22λm6){203λm2x3J3(x)5!!-[105-73λm2(6i1-5)]x4J4(x)7!!+910[7(60-48i1)+2λm2(8i2-7i1)]x5J5(x)9!!-9912[(35-112i1+80i2)+2λm2i2]x6J6(x)11!!-9·11·1314(7i1-20i2)x7J7(x)13!!-9·11·13·1516·i2x8J8(x)15!!}.
Substituting (26) and (30) into (6), the solution of the dimensionless luminal concentration profile is therefore(48)C(θ,x)=∑m=1∞∑n=0∞Fm(θ)Tmn(x)κn.
The perturbation solution is only extended up to second-order approximation; therefore, (48) reduces to(49)C(θ,x)=∑m=1∞Fm(θ)[Tm0(x)κ0+Tm1(x)κ1+Tm2(x)κ2].


## 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(50)u∂c∂z-DAB∂2c∂z2=DABr(∂c∂r+r∂2c∂r2)-v∂c∂r-VmcKm+c,
where Vm is the maximum rate of reaction and Km 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 Km*, respectively:(51)Km*=Kmc0,ϕ=VmRL2c0DAB.
Equation (50) then becomes(52)φPeuU∂C∂Z-φ2∂2C∂Z2=1R(∂C∂R+R∂2C∂R2)-PevV∂C∂R-ϕ2CKm*+C.
Assuming the first-order limit of the Michaelis-Menten equation (i.e., Km*≫C) 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 Fm(θ) and the coefficient of the perturbation solution B1m to account for the reaction rate:(53)C(θ,x)=∑m=1∞F̃m(θ)[Tm0(x)κ0+Tm1(x)κ1+Tm2(x)κ2],
where(54)F̃m(θ)=F0M(A(λm2Km*+ϕ2)2φ2Km*,12,θ),B̃1m=2λmF0M(A(λm2Km*+ϕ2)/2φ2Km*,1/2,θ0)×[J1(λm)J02(λm)+J12(λm)].
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 (km), the fraction retentate (f), the lumen-side inlet hydrostatic pressure (p0), the shell-side hydrostatic pressure (pS), 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.Table 2
Parameter values used to determine the concentration profile [17].
Model parameterSymbolUnitBasic measured valueMembrane hydraulic permeabilityakmm/Pas3.82× 10−11Fraction retentateaF0.70Membrane inner radiusRLm1.30× 10−4Annulus radiusr3m4.08× 10−4Effective membrane lengthLm5.7× 10−2Lumen-side entrance axial velocityu0ms−11.67× 10−3Permeation velocityv0ms−13.82× 10−7Number of fibresN150Lumen-side inlet hydrostatic pressureap0Pa106 325Shell-side hydrostatic pressureapSPa101 325Glucose diffusivityDABm2s1.0× 10−10solution densityaρkg m−3998solution viscosityaμPas9.7× 10−4Glucose inlet concentrationc0g dm−32.00Kinetic constantsVm/Kms−11.00aTypical operational values are used.In Figure3, 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 (R=0) 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 (R=1) 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 R 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].Figure 3
A comparison of the concentration profiles resulting from (53) with Heath and Belfort [17] assuming a first-order limit for substrate consumption (at R=0).Figure 4
Substrate concentration profiles assuming a first-order limit for substrate consumption at different radial positions.


## 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(50)u∂c∂z-DAB∂2c∂z2=DABr(∂c∂r+r∂2c∂r2)-v∂c∂r-VmcKm+c,
where Vm is the maximum rate of reaction and Km 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 Km*, respectively:(51)Km*=Kmc0,ϕ=VmRL2c0DAB.
Equation (50) then becomes(52)φPeuU∂C∂Z-φ2∂2C∂Z2=1R(∂C∂R+R∂2C∂R2)-PevV∂C∂R-ϕ2CKm*+C.
Assuming the first-order limit of the Michaelis-Menten equation (i.e., Km*≫C) 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 Fm(θ) and the coefficient of the perturbation solution B1m to account for the reaction rate:(53)C(θ,x)=∑m=1∞F̃m(θ)[Tm0(x)κ0+Tm1(x)κ1+Tm2(x)κ2],
where(54)F̃m(θ)=F0M(A(λm2Km*+ϕ2)2φ2Km*,12,θ),B̃1m=2λmF0M(A(λm2Km*+ϕ2)/2φ2Km*,1/2,θ0)×[J1(λm)J02(λm)+J12(λm)].
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 (km), the fraction retentate (f), the lumen-side inlet hydrostatic pressure (p0), the shell-side hydrostatic pressure (pS), 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.Table 2
Parameter values used to determine the concentration profile [17].
Model parameterSymbolUnitBasic measured valueMembrane hydraulic permeabilityakmm/Pas3.82× 10−11Fraction retentateaF0.70Membrane inner radiusRLm1.30× 10−4Annulus radiusr3m4.08× 10−4Effective membrane lengthLm5.7× 10−2Lumen-side entrance axial velocityu0ms−11.67× 10−3Permeation velocityv0ms−13.82× 10−7Number of fibresN150Lumen-side inlet hydrostatic pressureap0Pa106 325Shell-side hydrostatic pressureapSPa101 325Glucose diffusivityDABm2s1.0× 10−10solution densityaρkg m−3998solution viscosityaμPas9.7× 10−4Glucose inlet concentrationc0g dm−32.00Kinetic constantsVm/Kms−11.00aTypical operational values are used.In Figure3, 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 (R=0) 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 (R=1) 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 R 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].Figure 3
A comparison of the concentration profiles resulting from (53) with Heath and Belfort [17] assuming a first-order limit for substrate consumption (at R=0).Figure 4
Substrate concentration profiles assuming a first-order limit for substrate consumption at different radial positions.


## 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 (f); the membrane hydraulic permeability (κ); the axial and radial Peclet numbers (Peu,v); 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.



---

*Source: 738482-2010-08-24.xml*