Finite element model for nutrient distribution analysis of a hollow fiber membrane bioreactor G. U. Unnikrishnan 1 , V. U. Unnikrishnan 2 and J. N. Reddy 2, * ,† 2 Advanced

1 Department of Mechanical Engineering, Boston University, Boston, MA 02215, USA Computational Mechanics Laboratory, Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843, USA

SUMMARY Hollow fiber membrane bioreactors (HFMB) are extensively used for the development of tissue substitutes for bones and cartilages. In an HFMB, the nutrient transport is dependent on the material properties of the porous scaffold and fiber membrane and also on the fluid flow through the hollow fiber. The difficulty in obtaining real-time data along with the presence of large number of variables in experimental studies have lead to increased application of computational models for the performance analysis of bioreactors. A major difficulty in the computational analysis of HFMB is the modeling of the interactions at the fluid and porous scaffold interfaces, which has often been neglected or incorporated using specific boundary conditions. In this study, a new FEM is developed to analyze the fluid flow in the fluid-porous region with the interface coupled directly into the FEM. Distribution of nutrients in the bioreactor is studied by coupling mass transport equations to the fluid-porous finite element framework. The new model is implemented to study the influence of permeability, cell density, and flow rate on the nutrient concentration distribution in the HFMB. The developed computational framework is an ideal tool to study fluid flow through porous-open channels and can also be used for the design of bioreactors for optimal tissue growth. Copyright © 2011 John Wiley & Sons, Ltd. Received 18 February 2011; Revised 15 June 2011; Accepted 16 June 2011 KEY WORDS:

bioreactor; theory of mixtures; fluid-porous interface; penalty finite element method

1. INTRODUCTION Bioreactors simulate the physiological environment necessary for the development of tissue substitutes and play an important role in tissue engineering. In a bioreactor, the cells are attached to the porous scaffold that provides mechanical support and shape to the final tissue substitute [1–4]. The bioreactor maintains optimal growth environments by controlling the temperature, pH, and nourishment as in an in-vivo environment [5]. The major operational difficulty in the working of a bioreactor is in maintaining adequate nourishment and also removing waste materials from the scaffold [3, 6, 7]. Hollow fiber membrane bioreactor (HFMB), shown in Figure 1, which consists of a network of semi-permeable hollow fibers embedded in a porous scaffold, are widely used in the tissue engineering of bones and cartilages [8] because of the excellent nutrient transfer capabilities. In an HFMB, the transport and uptake of the nutrients by the cells in the scaffold are governed by a series of convection-diffusion-reaction mechanisms. The fluid flowing through the hollow fiber, material properties of the fiber and scaffold and also fluid-scaffold interactions affect the distribution of nutrient in the bioreactor [2, 3, 6, 9–12].

*Correspondence to: J. N. Reddy, J. N. Reddy, Advanced Computational Mechanics Laboratory, Department of Mechanical Engineering, Texas A&M University, 3123 TAMU College Station TX, USA. † E-mail: [email protected] Copyright © 2011 John Wiley & Sons, Ltd.

230

G. U. UNNIKRISHNAN, V. U. UNNIKRISHNAN AND J. N. REDDY scaffold matrix

hollow fiber membrane

Outlet

Inlet

direction of fluid flow

Figure 1. Schematic representation of hollow fiber membrane bioreactor (HFMB) of porous medium.

Fluid flow and nutrient transport analysis in HFMB are carried out by experiments, analytical solutions, and by numerical methods [7, 13, 14]. The presence of large number of variables and the difficulty in obtaining real time data makes experimental studies tedious. Numerical studies offer a viable alternative to study the complex fluid flow and scaffold interactions in the bioreactor. To simulate the uptake of solute in the bioreactor, the fluid flow and mass transport equations should be concurrently solved in the fluid and scaffold/fiber domain. Development and implementation of proper matching boundary conditions satisfying the continuity of mass, energy, and momentum at the fluid-porous interface have lead to models solving fluid and mass transport equations in the individual domains. In these two-domain models, additional boundary conditions are often explicitly defined in the numerical solution. In contrast to the two-domain model, a single-domain model treats the entire domain as a single continuum, and the same governing equations are used to model the fluid flow in the fluid and porous regions. The single-domain approach requires minimal alteration to the numerical modeling and is efficient for problems involving complex geometric interfaces. The major objective of this work is to develop an FEM for the simulation of fluid flow and nutrient transport in the HFMB using a single domain approach. A theory of mixtures framework, which assumes coexistence of different phases at a point in space [15, 16], is adopted to develop the singledomain FEM. In the FEM, the scaffold-fiber-fluid domain is treated as a single continuum with the interface coupled automatically into the governing differential equation. The scaffold and the fiber are treated as the solid phase, whereas the fluid culture medium flowing through the hollow fiber and the scaffold forms the fluid phase. The use of interfacial elements and requirement of any special treatment of the boundary conditions is avoided in the finite element formulation. The mass transfer in the bioreactor is then studied by incorporating convection-diffusion equations into the theory of mixtures framework. The nutrient considered in this study is oxygen as it is one of the critical factors in the growth of tissues [3, 7, 8]. Nutrient transfer is assumed to occur by both convection and diffusion in the scaffold, which is a deviation from previous works on nutrient transfer [3, 8]. Following this introduction, the governing equations of the theory of mixtures and the FEM are presented in “Section 2”. In “Section 3”, verification of the fluid-porous FEM is carried out. Also, analysis of nutrient distribution and fluid flow in the hollow fiber membrane bioreactor under different material properties are studied in this section. Conclusions of the study are presented in “Section 4”. 2. DEVELOPMENT OF FINITE ELEMENT MODEL In this section, the development of the theory of a mixture-based FEM for the analysis of fluid-porous domain analysis is presented. 2.1. Conservation of mass and momentum Conservation of mass of the different constituents in the continuum is given by Copyright © 2011 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Biomed. Engng. 2012; 28:229–238 DOI: 10.1002/cnm

FEM FOR NUTRIENT DISTRIBUTION ANALYSIS OF AN HFMB

231

@˛ C r. .˛ v˛ / D Γ˛ @t

(1)

where ˛ is the density, v˛ is the velocity and Γ˛ is the mass source/sink term for the ˛th constituent, which may depend on the volume fraction of the constituent. In the analysis of mixture containing n phases, the mass source/sink term should satisfy the constraint n X

Γ˛ D 0

(2)

˛D1

The balance of momentum of the individual phases can be represented as @ .˛ ˛ v˛ / C r. .˛ ˛ v˛ v˛ / D r.T˛ C ˛ ˛ b˛ C M˛ C Γ˛ vki @t

(3)

where T˛ is the stress tensor, b˛ is the body force, M˛ is the linear momentum because of interaction with other components inside the body, and Γ˛ vki is the momentum term because of the source of mass. Similar to the mass source/sink term constraint, the momentum interaction constraint should satisfy the condition n X

.M˛ C Γ˛ vki / D 0

(4)

˛D1

2.2. Constitutive equations of phases The constitutive expressions of the phases are obtained by simply generalizing the structure of the constitutive relations from a single constituent theory. In general, the constitutive equations for the stress Ts , Tf depend on the kinematical quantities of both the constituents, but in this work, we assume the principle of phase separation, where the stresses are assumed to depend only on the kinematic quantities associated with the individual components. The following simplifications are considered in this study: 1. The solid phase (scaffold and fiber) is a rigid porous medium with a constant porosity. 2. Fluid phase in the domain is viscous and incompressible. 3. Interaction term exists only because of the frictional term, and it is proportional to the difference in the velocity of phase components. The interaction is given by Mf D Cf vs vf (5) where Cf is the constant drag coefficient and vs D 0. The final governing equations for the fluid flow in the fluid and porous domains are given by r. pI C 2f Df C Cf vf D 0 r.vf D 0

(6) (7)

The interaction term Cf D 0 for the flow through the fluid domain, and in the rigid porous domain, the interaction term is given by Cf D

1 0

(8)

where 0 denotes the permeability constant. It is to be noted that unlike the Darcy or Brinkman equation, the interaction terms based on mixture theory is capable of accounting for complex interactions because of phase diffusion, virtual mass effect [15]. Copyright © 2011 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Biomed. Engng. 2012; 28:229–238 DOI: 10.1002/cnm

232

G. U. UNNIKRISHNAN, V. U. UNNIKRISHNAN AND J. N. REDDY

2.3. Mass transfer in porous system The mass transport of the solute.ˇ/, within the lumen and porous scaffold material is given by the following convective-diffusion equation r. vf C ˇ D ˇ rC ˇ D q ˇ (9) where, C ˇ (mol/m 3 / is the concentration of the solute in the medium, vf is the velocity of the fluid phase, q ˇ is the source term that considers the generation, consumption, or degradation of the solute mass, D ˇ is the effective diffusive coefficient in the porous medium. The diffusion coefficient is dependent on the molecular size, the solvent, and temperature [17, 18], and as the molecular size increases, the diffusion coefficient decreases. The source term in Equation (9) is represented by a simple Michaelis–Menton relationship for nutrient metabolism, and is given as [4, 19] V ˇ cell Cˇ qˇ D ˇ Km C C

(10)

ˇ (mol/m3 / is the concentration at which where V ˇ (mol/cell.s) is the maximum uptake rate, Km cell 3 the half-maximum rate occurs, and (cell/m / is the density of cell seeded in the scaffold. For ˇ ˇ Km