# A Finite-Difference Solution of Solute Transport through a Membrane Bioreactor

Authors: B. Godongwana; D. Solomons; M. S. Sheldon
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The current paper presents a theoretical analysis of the transport of solutes through a fixed-film membrane bioreactor (MBR), immobilised with an active biocatalyst. The dimensionless convection-diffusion equation with variable coefficients was solved analytically and numerically for concentration profiles of the solutes through the MBR. The analytical solution makes use of regular perturbation and accounts for radial convective flow as well as axial diffusion of the substrate species. The Michaelis-Menten (or Monod) rate equation was assumed for the sink term, and the perturbation was extended up to second-order. In the analytical solution only the first-order limit of the Michaelis-Menten equation was considered; hence the linearized equation was solved. In the numerical solution, however, this restriction was lifted. The solution of the nonlinear, elliptic, partial differential equation was based on an implicit finite-difference method (FDM). An upwind scheme was employed for numerical stability. The resulting algebraic equations were solved simultaneously using the multivariate Newton-Raphson iteration method. The solution allows for the evaluation of the effect on the concentration profiles of (i) the radial and axial convective velocity, (ii) the convective mass transfer rates, (iii) the reaction rates, (iv) the fraction retentate, and (v) the aspect ratio.


---

## Body

## 1. Introduction

Membrane bioreactors (MBRs) are finding increasing use in the production of primary and secondary metabolites such as amino acids, antibiotics, anticancer drugs, and tissue cells [1–3]. This technology is favoured by recent trends towards environmentally-friendly technologies, particularly because MBRs do not require additives, function at moderate operating conditions, and reduce by-product formation [1]. The efficiency of MBRs is dependent mainly on the transport of solutes through the bioreactor, and this is influenced by biochemical, geometric, and hydrodynamic parameters [2, 4]. This paper considers the numerical solution of the convection-diffusion equation for solute transport through a fixed-film MBR. This analysis is important for simulation of the performance (i.e., efficiency and effectiveness) of the bioreactor. The governing equation for mass transport of solutes through the bioreactor is the convection-diffusion equation with Monod kinetics [5] as follows:(1)u∂c∂z+v∂c∂r=DAB1r∂∂rr∂c∂r+∂2c∂z2-VMcKm+c,where c is the local substrate concentration, u and v are the axial and radial velocity components, respectively, DAB is the substrate diffusion coefficient, Vm is the maximum rate of reaction, and Km is the saturation (or Michaelis) constant. Equation (1) is made dimensionless by introducing the following variables:(2)U=uu0;V=vv0;C=cc0;ϕ=VMR12c0DAB,Km*=Kmc0;Z=zL;R=rR1;φ=R1L.Equation (1) then becomes(3)φPeuU∂C∂Z+PevV∂C∂R=1R∂∂RR∂C∂R+φ2∂2C∂Z2-ϕ2CKm*+C,where the axial and radial Peclet numbers Peu,v are, respectively, defined as(4)Peu=u0R1DAB;Pev=v0R1DAB.The velocity profiles, U and V, in (3) are solutions of the z and r-components of the Navier-Stokes equations, respectively [6]:(5)U=-141-R2dPdZ-ReFr,V=φu0v0R81-R22d2PdZ2,where P is the dimensionless hydrostatic pressure which is a function of the membrane hydraulic permeability κ and Re and Fr are the Reynolds and Froude numbers, respectively. When the membrane hydraulic permeability κ is much smaller than unity, (3) reduces to(6)U*∂C∂Z+Peuκβ2R-R3∂C∂R=1R∂∂RR∂C∂R+φ2∂2C∂Z2-ϕ2CKm*+C,where (7)U*=-4Peuκβ1f-1+Z1-R2,f≠1.The fraction retentate, f, is defined as the ratio of the outlet to the inlet axial velocity (f=0 for the dead-end mode and f~1 for the closed-shell mode), and β is the dimensionless transmembrane pressure. The corresponding boundary conditions are(8)B.C.1atZ=0∀RC=1,B.C.2atR=0∀Z∂C∂R=0,B.C.3atR=1∀Z∂C∂Z=2φPeu∂C∂R.Boundary condition 1 (B.C.1) corresponds to a uniform inlet substrate concentration; B.C.2 corresponds to cylindrical symmetry at the center of the membrane lumen; B.C.3 corresponds to continuity of the substrate flux at the lumen-matrix interface [7]. The solution of (6) is based on the following general assumptions: (i) the system is isothermal; (ii) the flow regime is laminar and fully developed; (iii) the fluid is Newtonian and homogenous and has constant physical and transport properties; (iv) the membrane hydraulic permeability is constant.A schematic of the MBR is shown in Figure1. The MBR 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. The nutrient solution is supplied by a peristaltic pump and permeates from the lumen-side to the shell-side of the MGR. The microorganism is immobilised on either the lumen-side or 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.Figure 1
A schematic diagram of the single capillary membrane bioreactor.


## 2. Analytical Models

The Graetz problem [8] is one of the oldest forced-convection problems describing the steady-temperature distribution and rate of heat transfer in tube flow. The evaluation of (6) for concentration profiles in a tubular reactor is mathematically analogous to the Graetz problem [5, 9]. In the original Graetz problem, however, there is no reaction (or source) term and axial diffusion, and radial convection is ignored. The assumption of negligible axial diffusion and radial convection is common in the majority of analytical models currently in use [7, 10–12]. Radial convective flows have been shown to significantly improve MBR efficiency [4, 13–15]. In the dead-end ultrafiltration mode, particularly, the assumption of negligible radial convective flow is not justifiable. Nagy [15] investigated the effect of radial convective flows (Pev) on the mass transfer rates of solutes through a biocatalytic membrane layer. Analytical solutions of (6) for the zero-order and first-order limits of the Monod equation were provided. This analysis, however, was restricted to the matrix/fiber region of the membrane and hence the radial velocity was assumed constant, and axial convective and diffusive flows were ignored. At axial Peclet numbers smaller than one (Peu≪1) large concentration gradients exist in the membrane lumen, and ignoring axial diffusion is also not justified [6, 7].The model proposed by Godongwana et al. [6] follows the approach suggested by Davis [9], that is, writing the solution of (6) in terms of known functions. The model accounts for radial convective flow and axial diffusion for the limiting case of first-order kinetics. In that model, (6) was solved by separation of variables and regular perturbation, resulting in the asymptotic expansion:(9)Cθ,x=∑m=1∞∑n=0NBmFmθTnxκn.Making use of the following changes of variables:(10)ξ=-2Peuκβφ21f-1+Z,θ=-φ24Peuκβξ2,x=λmR,where F(θ) in (9) is represented by the Kummer function [16] as follows:(11)Fmθ=M-λm2+ϕ2/Km*4Peuκβ,12,θ.The zero-order and first-order approximations of T(x) in (9) are, respectively,(12)T0x=J0x,T1x=σ1x2J2x3‼+σ2x3J3x5‼+σ3x4J4x7‼,where λm are the eigenvalues and Jn is the Bessel function of the first kind of order n [16]:(13)σ1=-3Pevβu0/v02λm2,σ2=-203λm2,σ3=354λm2.The eigenvalues are obtained from B.C.3 of (8) and are roots of the following equation:(14)φξ2κβλm2+ϕ2Km*M-λm2+ϕ2/Km*4Peuκβ+1,32,θ=2λmJ1λJ0λM-λm2+ϕ2/Km*4Peuκβ,12,θ.The first ten eigenvalues, at different axial positions, are listed in the appendix. The coefficient Bm is obtained by imposing the inlet condition B.C.1 of (8) and employing Lommel integrals to give [6](15)Bm=2λmM-λm2+ϕ2/Km*/4Peuκβ,1/2,θ0·J1λmJ02λm+J12λm.The assumption of first-order kinetics Km*≫C, used in the above analysis, allows for analytical evaluation of (6); however it limits the range of inlet substrate concentrations. The complete nonlinear form of (6) is not amenable to analytical evaluation and hence was solved using a finite-difference scheme described in Section 3.


## 3. Finite-Difference Scheme

A finite-difference representation of (6) is obtained by employing first-order upwind difference quotients for the derivatives on the LHS:(16)∂C∂Ri,j=Ci,j-Ci-1,jh+Oh,∂C∂Zi,j=Ci,j-Ci,j-1k+Ok.And second-order central differences for the derivatives on the RHS of (6) are(17)∂2C∂R2i,j=Ci+1,j-2Ci,j+Ci-1,jh2+Oh2,∂2C∂Z2i,j=Ci,j+1-2Ci,j+Ci,j-1k2+Ok2,resulting in the general difference equation as follows:(18)α1Ci-1,j+α2+α3Km*+Ci.jCi,j+α4Ci+1,j+α5Ci,j-1+α6Ci,j+1=0,where(19)α1=-hPeuκβ2R-R3+1-hR,α2=h2kU*+hPeuκβ2R-R3+2φ2h2k2+2-hR,α3=h2ϕ2,α4=-1,α5=-h2U*k+φ2k2,α6=-φ2h2k2.The corresponding difference quotients of the boundary conditions in (8) are(20)B.C.1atj=1∀iCi,j=e-PevVR,B.C.2ati=1∀jCi,j-Ci-1,jh=0,B.C.3ati=m∀jCi+1,j-Ci,jh=φPeu2Ci,j-Ci,j-1k,B.C.4atj=n∀iCi,j-Ci,j-1k=0.Boundary condition 1 is an approximation of the inlet condition in (8) and is obtained by assuming the net diffusive and convective flux in the radial direction are negligible at the membrane entrance. Boundary condition 4 is applicable in the dead-end mode (one end of the lumen-side is closed). The solution domain is a regular 2-dimensional grid and is subdivided into m-intervals (of size h) in the r-dimension and n-intervals (of size k) in the z-dimension. The difference equation (18), including the boundary conditions of (20), is solved by making use of the multivariate Newton-Rapshon iteration scheme:(21)Cn+1=Cn-FCnJn,where F(C) is the residual equation (18) and J is the tridiagonal Jacobian matrix as follows:(22)J=DE0BDE0BD⋯000000000⋮⋱⋮000000000⋯DE0BDE0BD.The matrix elements B, D, and E are, respectively,(23)B=α4+α5,D=α2+α31Km*+Ci.j+Ci.jKm*+Ci.j2,E=α1+α6.The Newton-Raphson iteration scheme was implemented on MATLAB R2014a and the procedure is shown in Figure 2. The algorithm begins with an initial guess of the solute concentration at each grid point; an initial guess of zero was used. Based on this guess, the residual column vector and the Jacobian matrix can be evaluated. The magnitude (Euclidean norm) of the quotient of the residual vector and Jacobian matrix, dC, is evaluated. The iteration is repeated with new solute concentration guess values until the Euclidean norm is less than the prescribed tolerance.Figure 2
The Newton-Raphson algorithm for solving (18).


## 4. Results

### 4.1. Numerical Solution

The implicit finite-difference scheme was shown to be unconditionally stable for the different values ofh and k listed in Table 1. The results are shown in Figure 3 for the parameter values listed in Table 2. Figure 3 illustrates the effect of the fraction retentate f on the solute concentration profiles. In the dead-end mode (f=0) there is increased radial convective flow as shown by the streamlines in Figure 3(a). This increased radial flow allows for more solute contact with the biofilm and hence improved conversion, resulting in higher MBR efficiency. In this mode however the solute is limited to only the entrance half of the MBR as shown in Figure 3(a). Increasing the fraction retentate to f=0.8 allows for a uniform distribution of the solute (Figure 3(b)); however radial convective flow is significantly reduced. This result implies that an optimum f value should be sought for enhanced MBR efficiency. The developed finite-difference scheme also allows for the evaluation of the effect on the concentration profiles of the radial and axial convective velocity, the convective mass transfer rates, the reaction rates, and the aspect ratio. The sensitivity analysis of these parameters however has been omitted in the current paper.Table 1
Computation times for different sizes ofh and k.
Simulation
Number of radial nodes (m)
Number of axial nodes (n)
Total number of solution nodes
Total iteration time (s)
1
32
32
1,024
13.302
2
64
64
4,096
55.130
3
128
128
16,384
505.108Table 2
Parameter values used to determine the concentration profile [6].
Model parameter
Symbol
Unit
Basic measured value
Membrane hydraulic permeability
k
m
m/Pas
3.82× 10−11
Membrane inner radius
R
1
m
1.30× 10−4
Effective membrane length
L
m
5.7× 10−2
Lumen-side entrance axial velocity
u
0
m s−1
1.67× 10−3
Permeation velocity
v
0
m s−1
1.91× 10−7
Lumen-side inlet hydrostatic pressure
p
0
Pa
106 325
Shell-side hydrostatic pressure
p
S
Pa
101 325
Glucose diffusivity
D
A
B
m2 s
1.0× 10−10
Solution density
ρ
kg m−3
998.0
Solution viscosity
μ
Pas
9.7× 10−4
Glucose inlet concentration
c
0
g dm−3
2.00
Kinetic constants
K
m
/
V
m
s−1
0.10Figure 3
Solute concentration profiles form=n=64 when (a) f=0 and (b) f=0.8 (for the parameter values listed in Table 2).
(a)
(b)Figure4 is a plot of the solute concentration profile, from the FDM scheme, for different dimensionless saturation constants when f=0.8. At a low Km* value of 0.009 the dimensionless outlet solute concentration is 0.472, whereas increasing Km* to 0.9 reduces the outlet concentration to only 0.853. This result is expected since a lower Km* value is consistent with a higher biofilm/enzyme affinity for the substrate [17]. This result also illustrates the significance of the appropriate choice of substrate on bioreactor design.Figure 4
Solute concentration profiles from the finite-difference method (FDM) for different saturation constantsKm, when f=0.8.


### 4.2. Comparison with Analytical Solution

In Figure5(a) the FDM scheme is compared with the analytical model presented in Section 2 for the open-shell mode (f=0.8). The analytical model predicts a linear decrease in the solute concentration inside the membrane lumen to 47% of the original concentration. This result is consistent, qualitatively, with the result of Heath and Belfort [7] for the parameter values listed in Table 2. The FDM scheme predicts the same outlet concentration; however the decrease is gradually close to the entrance and rises with increasing length. The discrepancies between the two profiles arise from the assumption of first-order kinetics, assumed in developing the analytical solution. These two profiles suggest that the open-shell mode is suitable for microbial growth since the substrate is not depleted inside the lumen.Figure 5
A comparison of the analytical versus the FDM solution for solute concentration profiles: (a)f=0.8 and (b) f=0 (for the parameter values listed in Table 2).
(a)
(b)The rapid decline in the solute concentration in Figure5(b) is due to increased radial convective flow in the dead-end mode. This results in nonuniform microbial growth/tapering as observed by Godongwana et al. [18] for the bacteriumStreptomyces coelicolor on a ceramic membrane. This phenomenon can be reduced by either increasing the solute flow rate or increasing the fraction retentate f. The numerical scheme matches the analytical model approximately on a small interval close to the origin. The divergence again is attributed to the assumption of first-order kinetics.


## 4.1. Numerical Solution

The implicit finite-difference scheme was shown to be unconditionally stable for the different values ofh and k listed in Table 1. The results are shown in Figure 3 for the parameter values listed in Table 2. Figure 3 illustrates the effect of the fraction retentate f on the solute concentration profiles. In the dead-end mode (f=0) there is increased radial convective flow as shown by the streamlines in Figure 3(a). This increased radial flow allows for more solute contact with the biofilm and hence improved conversion, resulting in higher MBR efficiency. In this mode however the solute is limited to only the entrance half of the MBR as shown in Figure 3(a). Increasing the fraction retentate to f=0.8 allows for a uniform distribution of the solute (Figure 3(b)); however radial convective flow is significantly reduced. This result implies that an optimum f value should be sought for enhanced MBR efficiency. The developed finite-difference scheme also allows for the evaluation of the effect on the concentration profiles of the radial and axial convective velocity, the convective mass transfer rates, the reaction rates, and the aspect ratio. The sensitivity analysis of these parameters however has been omitted in the current paper.Table 1
Computation times for different sizes ofh and k.
Simulation
Number of radial nodes (m)
Number of axial nodes (n)
Total number of solution nodes
Total iteration time (s)
1
32
32
1,024
13.302
2
64
64
4,096
55.130
3
128
128
16,384
505.108Table 2
Parameter values used to determine the concentration profile [6].
Model parameter
Symbol
Unit
Basic measured value
Membrane hydraulic permeability
k
m
m/Pas
3.82× 10−11
Membrane inner radius
R
1
m
1.30× 10−4
Effective membrane length
L
m
5.7× 10−2
Lumen-side entrance axial velocity
u
0
m s−1
1.67× 10−3
Permeation velocity
v
0
m s−1
1.91× 10−7
Lumen-side inlet hydrostatic pressure
p
0
Pa
106 325
Shell-side hydrostatic pressure
p
S
Pa
101 325
Glucose diffusivity
D
A
B
m2 s
1.0× 10−10
Solution density
ρ
kg m−3
998.0
Solution viscosity
μ
Pas
9.7× 10−4
Glucose inlet concentration
c
0
g dm−3
2.00
Kinetic constants
K
m
/
V
m
s−1
0.10Figure 3
Solute concentration profiles form=n=64 when (a) f=0 and (b) f=0.8 (for the parameter values listed in Table 2).
(a)
(b)Figure4 is a plot of the solute concentration profile, from the FDM scheme, for different dimensionless saturation constants when f=0.8. At a low Km* value of 0.009 the dimensionless outlet solute concentration is 0.472, whereas increasing Km* to 0.9 reduces the outlet concentration to only 0.853. This result is expected since a lower Km* value is consistent with a higher biofilm/enzyme affinity for the substrate [17]. This result also illustrates the significance of the appropriate choice of substrate on bioreactor design.Figure 4
Solute concentration profiles from the finite-difference method (FDM) for different saturation constantsKm, when f=0.8.


## 4.2. Comparison with Analytical Solution

In Figure5(a) the FDM scheme is compared with the analytical model presented in Section 2 for the open-shell mode (f=0.8). The analytical model predicts a linear decrease in the solute concentration inside the membrane lumen to 47% of the original concentration. This result is consistent, qualitatively, with the result of Heath and Belfort [7] for the parameter values listed in Table 2. The FDM scheme predicts the same outlet concentration; however the decrease is gradually close to the entrance and rises with increasing length. The discrepancies between the two profiles arise from the assumption of first-order kinetics, assumed in developing the analytical solution. These two profiles suggest that the open-shell mode is suitable for microbial growth since the substrate is not depleted inside the lumen.Figure 5
A comparison of the analytical versus the FDM solution for solute concentration profiles: (a)f=0.8 and (b) f=0 (for the parameter values listed in Table 2).
(a)
(b)The rapid decline in the solute concentration in Figure5(b) is due to increased radial convective flow in the dead-end mode. This results in nonuniform microbial growth/tapering as observed by Godongwana et al. [18] for the bacteriumStreptomyces coelicolor on a ceramic membrane. This phenomenon can be reduced by either increasing the solute flow rate or increasing the fraction retentate f. The numerical scheme matches the analytical model approximately on a small interval close to the origin. The divergence again is attributed to the assumption of first-order kinetics.


## 5. Conclusion

A numerical solution of the dimensionless convection-diffusion equation, with nonlinear kinetics, was developed. The numerical scheme was performed using the Newton-Raphson method and was shown to be unconditionally stable for different step-sizes (h and k). The analysis provides for evaluation of concentration profiles of solutes through a membrane bioreactor. The numerical solution was compared to a regular perturbation solution for two modes of operation, that is, the dead-end mode and open-shell mode. In the dead-end mode the numerical results closely matched the perturbation solution. The assumption of linear kinetics, commonly used in literature models, was shown to result in inaccuracies in the open-shell mode. The numerical solution allows for the evaluation of the influence of the general operating parameters of a MBR on the concentration profiles. The fraction retentate (f) was shown to be an important optimisation parameter for improved MBR efficiency.



---

*Source: 810843-2015-03-25.xml*