- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modelling of fluid flow and protein transport in hollow-fibre...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modelling of fluid flow and protein transport in hollow-fibre bioreactors 1994
pdf
Page Metadata
Item Metadata
Title | Modelling of fluid flow and protein transport in hollow-fibre bioreactors |
Creator |
Labecki, Marek |
Date Created | 2009-03-03 |
Date Issued | 2009-03-03 |
Date | 1994 |
Description | A mathematical model (the Porous Medium Model, PMM) was developed to predict the
fluid flow and solute transport in hollow-fibre devices, with a particular emphasis on hollowfibre
bioreactors (HFBRs). In the PMM, both the extracapillary space (ECS) and the lumen
side are treated as interpenetrating porous regions with a continuous source or sink of fluid.
The hydrodynamic equations of the PMM are based on Darcy's law and continuity
considerations while the transport of the ECS protein is described by the time-dependent
convective-diffusion equation. Compared to the earlier Krogh Cylinder Model (KCM), in
which the fluid flow and protein transport are assumed to be the same for each fibre, the
PMM represents an improved approach in which the spatial domain corresponds to the real
dimensions of the hollow-fibre module. Thus, it can be applied to operating conditions where
macroscopic radial pressure and concentration gradients exist, such as in open-shell
operations. It was demonstrated that, in the absence of radial gradients, the PMM becomes
mathematically equivalent to the one-dimensional KCM. The PMM also takes into account
the osmotic pressure dependence on the ECS protein concentration, which causes a coupling
of the hydrodynamic and protein transport equations.
The Porous Medium Model was tested by applying it to one- and two-dimensional
closed-shell operations. Both confirmed that a significant polarization of the ECS protein
occurs in the direction of the existing pressure gradients under dominant convective transport
conditions. The downstream polarization of protein affects HFBR hydrodynamics by virtually
shutting down the flow in a significant portion of the ECS due to locally high osmotic
pressures. It can also facilitate harvesting of the product protein by increasing its
concentration near the downstream ECS port. Modelling studies of the hydrodynamics of hollow-fibre devices in the partial and full
filtration modes of operation were carried out for a wide range of membrane permeabilities
(10"14 |
Extent | 5550436 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-03-03 |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0087507 |
Degree |
Master of Applied Science - MASc |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1994-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/5393 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0087507/source |
Download
- Media
- ubc_1994-0479.pdf [ 5.29MB ]
- Metadata
- JSON: 1.0087507.json
- JSON-LD: 1.0087507+ld.json
- RDF/XML (Pretty): 1.0087507.xml
- RDF/JSON: 1.0087507+rdf.json
- Turtle: 1.0087507+rdf-turtle.txt
- N-Triples: 1.0087507+rdf-ntriples.txt
- Citation
- 1.0087507.ris
Full Text
MODELLING OF FLUID FLOW AND PROTEIN TRANSPORT IN HOLLOW-FIBRE BIOREACTORS by MAREKLABSCKI M.Sc, Technical University of Szczecin, Poland, 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Chemical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA © MarekLabecki, May 1994 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Chemical Engineering University of British Columbia 2216 Main Mall Vancouver, B.C. V6T 1Z4 Canada Date: May 1994. ii ABSTRACT A mathematical model (the Porous Medium Model, PMM) was developed to predict the fluid flow and solute transport in hollow-fibre devices, with a particular emphasis on hollow- fibre bioreactors (HFBRs). In the PMM, both the extracapillary space (ECS) and the lumen side are treated as interpenetrating porous regions with a continuous source or sink of fluid. The hydrodynamic equations of the PMM are based on Darcy's law and continuity considerations while the transport of the ECS protein is described by the time-dependent convective-diffusion equation. Compared to the earlier Krogh Cylinder Model (KCM), in which the fluid flow and protein transport are assumed to be the same for each fibre, the PMM represents an improved approach in which the spatial domain corresponds to the real dimensions of the hollow-fibre module. Thus, it can be applied to operating conditions where macroscopic radial pressure and concentration gradients exist, such as in open-shell operations. It was demonstrated that, in the absence of radial gradients, the PMM becomes mathematically equivalent to the one-dimensional KCM. The PMM also takes into account the osmotic pressure dependence on the ECS protein concentration, which causes a coupling of the hydrodynamic and protein transport equations. The Porous Medium Model was tested by applying it to one- and two-dimensional closed-shell operations. Both confirmed that a significant polarization of the ECS protein occurs in the direction of the existing pressure gradients under dominant convective transport conditions. The downstream polarization of protein affects HFBR hydrodynamics by virtually shutting down the flow in a significant portion of the ECS due to locally high osmotic pressures. It can also facilitate harvesting of the product protein by increasing its concentration near the downstream ECS port. iii Modelling studies of the hydrodynamics of hollow-fibre devices in the partial and full filtration modes of operation were carried out for a wide range of membrane permeabilities (10"1 4<LP< 10-7 m). It was demonstrated using the PMM that, for membranes with permeabilities below about 10"13 m, practically all of the pressure drop between the inlet lumen and outlet ECS ports is due to the hydraulic resistance of the membrane. If the Lp value is increased above approximately 10'12 m, this assumption, commonly made in order to experimentally determine membrane permeabilities, begins to break down. Also, for membrane permeabilities exceeding this value, the ECS and lumen flow rates predicted by the PMM and KCM for the partial filtration mode become significantly different. Modelling of the inoculation phase of HFBR operation is used as another example application of the Porous Medium Model. PMM simulations of the inoculation phase showed that, in the case of a Gambro HFBR with a membrane permeability of the order of 10"15 m, the protein concentration distribution at the end of the inoculation period is very non-uniform and most of the shell side remains free of protein. Using a lower-concentration inoculum solution partially alleviates this problem. Alternatively, a relaxation phase with all ports closed can be applied after inoculation to help homogenize the contents of the ECS by diffusion and osmotically-driven convection. However, this process may be fairly time-consuming and may pose the risk of cell starvation due to oxygen limitations. It is suggested that introduction of the inoculum through both ECS ports simultaneously or periodic changes of the flow direction may be more efficient ways of carrying out the inoculation process. The cell-packed conditions, which exist in the ECS during the production and harvesting phases of HFBR operation, can significantly decrease the ECS hydraulic conductivity and, to a lesser extent, the effective protein diffusivity due to a decrease in the ECS porosity. The ECS permeability value affects the magnitude of convective transport in the shell side and hence the rate of protein removal from the ECS and the product concentration in the harvested solution, thereby influencing the overall efficiency of the process. High-cell-density conditions in the IV ECS might not allow achievement of high product removal rates and product harvest concentrations. Two modes of harvesting, the closed-lumen mode (with only the two ECS ports open) and the standard mode (with only the downstream ECS port and both lumen ports open), were compared and showed no significant differences in their efficiencies. It was found that the downstream polarization of the ECS protein prior to harvesting can considerably improve the efficiency of this process. V TABLE OF CONTENTS ABSTRACT ii TABLE OF CONTENTS v LIST OF TABLES vii LIST OF FIGURES viii ACKNOWLEDGEMENTS xi 1. INTRODUCTION 1 2. PREVIOUS MODELLING WORK 2.1. Introduction 6 2.2. General assumptions 6 2.3. Major approaches to modelling 9 3. DEVELOPMENT OF THE POROUS MEDIUM MODEL 3.1. Introduction 15 3.2. Hydrodynamics 17 3.3. Protein transport 19 3.4. Modelling of the ECS and lumen hydraulic conductivities and protein effective difHisivities 23 3.5. Boundary conditions 30 4. NUMERICAL TECHNIQUES 4.1. Solution of the pressure equations 36 4.2. Solution of the convective-diffusion equation 40 4.3. General computational algorithm 41 5. TESTING AND APPLICATION OF THE POROUS MEDIUM MODEL 5.1. Closed-shell operation: one-dimensional case 47 5.2. Closed-shell operation: two-dimensional case with inlet and outlet radial lumen pressure gradients 54 5.3. Membrane permeability determination 57 5.4. Filtration hydrodynamics: comparison of the Porous Medium Model with the Krogh Cylinder Model 64 5.5. Inoculation and relaxation 70 5.6. Harvesting 81 vi 6. CONCLUSIONS AND FUTURE WORK 96 NOMENCLATURE 100 REFERENCES 105 Appendix A: FORMULATIONS OF HFBR MODELLING EQUATIONS BASED ON THE KROGH CYLINDER APPROXIMATION I l l Appendix B: KROGH CYLINDER EQUATIONS WITH PRESSURE BOUNDARY CONDITIONS IN THE CLOSED-SHELL, PARTIAL, AND FULL FILTRATION MODES 116 Appendix C: PATANKAR'S POWER-LAW SCHEME 119 Appendix D: SOURCE CODE IN FORTRAN 121 Vll LIST OF TABLES 1.1: Some applications of hollow-fibre devices 1 2.1: Parameters affecting the flow in hollow-fibre bioreactors 7 3.1: Theoretical expressions for the dimensionless Darcy constant k/R^ as a function of the solid fraction <p for a flow through an array of parallel cylinders (RM is the cylinder radius); k values are calculated for <p=0.5 and RM =10_ 4m. a. Flow parallel to the cylinders 26 b. Flow perpendicular to the cylinders 27 3.2: Boundary conditions common to all flow configurations 32 3.3: Boundary conditions specific for different flow configurations 34 5.1: Parameters used in the comparative study of the PMM and KCM (Amicon HFBR) 50 5.2: Summary of the steady-state runs with different hypothetical diffusivities for the PMM and KCM (Amicon HFBR) 53 5.3: Important parameters used in the membrane permeability determination study 59 5.4: Experimentally determined volumetric flows, Q, and pressure drops, AP, for a Gambro HFBR (Koska, 1993b) 61 5.5: Summary of parameters used in the inoculation study (Gambro HFBR) 71 5.6: The effect of osmotic pressure in the inoculation tests 72 5.7: Summary of parameters used in the harvesting study (Gambro HFBR) 83 5.8: Pressure drops corresponding to the initial ECS flow rate of 1.0000 cm3/min and (standard mode only) the initial lumen flow rate of 600.00 cm3/min 86 Vlll LIST OF FIGURES Schematic of a hollow fibre device (not to scale) 2 Krogh cylinder approximation: a) fibre arrangement throughout the reactor cross-section, b) longitudinal section of a single unit with major flow paths 11 Ring-shaped representative elementary volume (REV) with thickness Ar and length Ax 16 The REV cross-section diagram for the convective-diffusion equation 20 Diagram of the spatial domain boundaries in the Porous Medium Model 31 Cluster of control volumes contributing to the finite difference equations for a point (i,j) in the interior of the domain. Note: the lumen pressure equation has contributions from three cells only (constant j) 37 Block diagram of the general computational algorithm used in this study 42 The form of the pressure iteration block used in this study 43 Examples of the pressure iteration block algorithm (not used in this study): a) with the lumen pressure lagged behind the ECS pressure, b) with the lumen and ECS pressures iterated until both converge 44 General algorithm with fully coupled concentration and pressure equations 45 Radially-averaged protein concentration in the ECS as a function of axial position and time 51 Steady-state protein concentration profiles for different hypothetical diffusivities 52 Steady-state ECS protein concentration in the absence (a,c) and presence (b,d) of radial pressure gradients; (a,b) CQ = 10 kg/m3, (c,d) c0 = 20 kg/m3 56 Flow configuration for Lp determination 58 Volumetric flow as a function of AP: comparison of the Porous Medium Model with experiment 60 IX 5.6: The apparent (calculated) membrane permeability, Lp>app, versus the actual (input) membrane permeability, Lp. prediction of the Porous Medium Model 62 5.7: The ratio o£Lp>app/Lp versus the actual membrane permeability, Lp: prediction of the Porous Medium Model 63 5.8: Flow diagram for the partial filtration mode (jQUout >0) and full filtration mode {Qu „,=(>) 64 5.9: The lumen inlet volumetric flow as a function of membrane permeability in the partial and full filtration modes: comparison of the two models 66 5.10: The ECS and lumen outlet flow rates, Qs,out and Quout, as functions of membrane permeability in the partial filtration mode: comparison of the two models ...68 5.11: Flow diagram for inoculation phase of HFBR operation 70 5.12: ECS concentration field after a) 20 min, b) 40 min, c) 60 min inoculation, cin = 5 kg/m3 74 5.13: ECS concentration field after a) 2 min, b) 4 min, c) 6 min inoculation, c,„ = 50kg/m3 75 5.14: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field («£-108 m/s) after 1 h relaxation following 60 min of inoculation with cin = 5 kg/m3 77 5.15: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field («i-108 m/s) after 20 h relaxation following 60 min of inoculation with c,„ = 5 kg/m3 78 5.16: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field (uL-l0* m/s) after 1 h relaxation following 6 min of inoculation with cin = 50 kg/m3 79 5.17: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field (i/flO8 m/s) after 20 h relaxation following 6 min of inoculation with c/n = 50 kg/m3 80 5.18: Some of the possible modes of harvesting: a) standard, b) closed-lumen, c) with all ports open and equal pressures at the ECS and lumen inlets, d) with closed lumen inlet 82 X 5.19: Steady-state ECS protein concentration as a function of axial position for different ECS porosities 84 5.20: The ECS outlet flow as a function of the lumen flow for different ECS porosities in the standard ultrafiltration mode with Ps,dn — PL,N = 1.0 atm 85 5.21: Fraction of protein removed from the HFBR as a function of harvesting time at different ECS porosities (uniform initial concentration field) 87 5.22: The ECS outlet concentration and the concentration in the harvesting reservoir as functions of time (uniform initial concentration field, s'ECS = 26%) 88 5.23: ECS protein concentration field (kg/m3) in the closed-lumen harvesting with eEcs = 5%: a) after 5 min, b) after 30 min, c) after 60 min ,...89 5.24: ECS concentration field (kg/m3) after 2 h of harvesting (uniform initial field), a) e'ECS = 26%, closed-lumen, b) e'ECS = 5%, closed-lumen, c) s'ECS = 5%, standard 90 5.25: The ECS outlet concentration and the concentration in the harvesting reservoir as functions of time (both harvesting modes, uniform initial concentration field, 4cs = 5%) 92 5.26: Fraction of protein removed from the HFBR as a function of the total outflow from the ECS at different ECS porosities (both harvesting modes, polarized initial concentration field) 93 5.27: The ECS outlet concentration and the concentration in the harvesting reservoir as functions of the total outflow from the ECS at different ECS porosities (both harvesting modes, polarized initial concentration field) 94 5.28: a) ECS concentrations (kg/m3), b) ECS velocity vectors (in the central part of the ECS) and streamlines (near the port manifolds, where the magnitude of the flow is much larger than in the central part of the ECS), c) lumen velocities («z,-105 m/s) after 10 min of closed-lumen harvesting (s'ECS =5%, polarized initial concn. field) 95 xi ACKNOWLEDGEMENTS I wish to express my sincere gratitude to my supervisors, Drs. Bruce D. Bowen and James M. Piret, for their direction, support, and friendly assistance throughout this project. My grateful appreciation is due to Dr. B. D. Bowen for his extreme patience and kindly accommodating me within his work space whenever I needed it. I would also like to thank the other members of my committee, Dr. Joel Bert and Dr. Charles Haynes, who have agreed to devote their time and provide me with their expertise regarding this work. Many thanks are due to my landlady, Donelda, for creating a real home environment for me and an always excellent cuisine. Finally, I am grateful to my friends, Lu, Dale, Linda, Igor, Mohandes, Eugene, Leif, for their joyful appearances and sharing their vital energy. May the many others, mostly overseas, forgive me for not mentioning them by names. I feel particularly indebted to my mother, without whom everything would be completely different. She is the special person this work is dedicated to. Your worst enemy cannot harm you As much as your own thoughts, unguarded. But once mastered, No one can help you as much, Not even your father or your mother. The Buddha Chapter 1: Introduction 1 Chapter 1: INTRODUCTION A typical hollow-fibre module consists of a bundle of semi-permeable polymeric capillaries sealed inside a tubular cartridge. The device has found numerous applications in various fields, some of which are presented in Table 1.1. Besides those listed, other possibilities exist for its use in the food and fermentation industries, tanning and textile industries, in waste-water treatment, and other traditional fields where filtration or reverse osmosis is applied (Drioli, 1980; Michaels, 1980). Using the hollow-fibre module as an immobilized enzyme bioreactor was first proposed by Rony (1971), while Knazek et al. (1972) first reported using a hollow-fibre bioreactor for mammalian cell culture. Whole-cell Table 1.1: Some applications of hollow-fibre devices. Application Reference artificial kidney (hemodialysis) artificial pancreas hemofilters and hemodiafilters liver assist device hormone production monoclonal antibody production purification of biological macromolecules ultrapure water production water desalination Mahon (1960) Colton et al. (1980) Gohl & Konstantin (1986) Wolf (1980) Knazek et al. (1972) Piret&Cooney (1990a) Michaels (1980) Michaels (1980) Breslau et al. (1980), Hermans (1978) Chapter 1: Introduction 2 immobilization has important advantages over enzyme immobilization including elimination of the enzyme purification step and the ability of whole cells to catalyze multi-step reactions (Webster & Shuler, 1978). Protection of cells from shear stresses, high cell densities, and increased product concentrations are some advantages of immobilized cell cultures compared with suspension cultures (Piret & Cooney, 1990b). Hollow-fibre systems offer a particularly high surface-to-volume ratio, and thus high throughput capacity and high productivity. In the case of artificial organs, foreign cells can be immunoisolated (Kelsey et al., 1990). On the other hand, difficult sampling, nutrient and metabolite gradients, and scale-up limitations are some of the potential problems associated with hollow-fibre devices (Piret & Cooney, 1991; Piret et al., 1991). Figure 1.1: Schematic of a hollow fibre device (not to scale). Figure 1.1 shows a schematic diagram of a hollow-fibre bioreactor (HFBR). In a typical configuration, it contains cells packed to high densities in the extracapillary space (ECS) and thus physically separated from the major flow that enters and exits the reactor and passes Chapter 1: Introduction 3 through the fibre lumina. Low-molecular-weight nutrients (e.g., dissolved oxygen, glucose, etc.) permeate through the membrane into the ECS, while cell metabolic wastes are continuously removed from the ECS to the lumen flow. Macromolecular product proteins usually remain in the ECS but, depending on the size and shape of the molecule and on the membrane properties, might also migrate to the lumen. This flow configuration is known as the closed-shell mode (since both ECS ports are closed) and is commonly used in the production phase. Other variations of the closed-shell configuration include periodic alternation of the flow direction (Piret & Cooney, 1990a) or applying pulsatile flow at the inlet (Kim & Chang, 1983). Open-shell operation occurs during the inoculation phase, in which the cells are introduced into the shell side through the upstream ECS port, or the harvesting phase, in which the product is collected from the downstream ECS port. In the prediction of cell and product distributions in HFBRs, one must account for at least some and possibly all of the following phenomena: (i) diffusion, (ii) convection, (iii) osmosis, (iv) gravity, (v) adsorption, and (vi) metabolic reactions. Adsorption of cells or protein molecules can lead to membrane fouling, while osmotically active species, if present at sufficiently high concentrations, will influence transmembrane flows. Distributions of low- molecular species, for which the membrane is permeable, are usually affected by diffusion and metabolic reactions only (Piret & Cooney, 1991; Piret et al., 1991). However, for proteins, the influence of gravity and convection often cannot be neglected, a conclusion made by Piret and Cooney based on their evaluation of the Grashof number (ratio of buoyant forces to viscous resistance) and Peclet number (ratio of convective to diffusive transport) for typical process conditions. This has been confirmed by experimental observation of downstream polarization and sedimentation of cells and proteins (Piret & Cooney, 1990a). Convectively induced downstream polarization of protein may cause problems associated with ineffective use of the reactor space during the start-up phase of HFBR operation since the cell distribution has been shown to follow the growth factor distribution (Piret & Cooney, 1990a). Chapter 1: Introduction 4 On the other hand, it can be advantageous to polarize the protein in the downstream part of the reactor prior to harvesting in order to obtain a more concentrated product solution. Some authors (e.g., Pillarella & Zydney, 1990; Salmon et al., 1988) have claimed that increased ECS convective flow should improve the productivity of HFBRs. Thus, the analysis of HFBR convective fluid flow is an important step in the modelling, design and scale-up of these bioreactors. The mathematical models describing hydrodynamics as well as convective and diffusive protein transport in HFBRs have so far been based on the analysis of a single fibre unit assumed to be representative of the whole reactor (Apelblat et al., 1974; Kelsey et al., 1990; Pillarella & Zydney, 1990; Taylor et al., 1994). In the closed-shell mode, the pressure gradients between fibres are often small, and this assumption is fairly reasonable. However, it cannot lead to a realistic description of HFBRs in cases of (1) significant radial pressure variations in the lumen manifolds, (2) significant concentration gradients between fibres, caused by gravity or by protein entrapment in the ECS manifolds, or (3) for open-shell operations such as inoculation or harvesting (since the ECS manifolds are located at the reactor circumference). This work presents the development and gives example applications of a new model, in which the bundle of densely packed hollow fibres is treated as a porous medium and the spatial domain is determined by the real cartridge dimensions. This not only makes it possible to handle the open-shell as well as the closed-shell flow configurations but also allows for relatively straightforward extensions of the model to include, for instance, the presence of ECS or lumen manifolds or the influence of gravity on the hydrodynamics and protein distribution. Chapter 2 of this thesis presents a review of previous models of fluid flow and protein transport in HFBRs. The assumptions and mathematical development of the proposed Porous Medium Model (PMM) are given in Chapter 3, while Chapter 4 briefly describes the numerical techniques used. Section 1 of Chapter 5 presents examples of model verification Chapter 1: Introduction 5 through the comparison of some one-dimensional solutions with those obtained from single- fibre models. Section 5.2 presents solutions obtained using the PMM with imposed radial lumen pressure gradients at the HFBR inlet and outlet. Finally, predictions of the PMM with respect to the determination of membrane permeability, hydrodynamics in the filtration mode, as well as inoculation and harvesting operations are presented in subsequent sections of Chapter 5. The last chapter concludes the thesis and outlines possible future extensions to this model. Chapter 2: Previous modelling work 6 Chapter 2: PREVIOUS MODELLING WORK 2.1. Introduction Most mathematical models describing HFBRs start with a hydrodynamic analysis and then superimpose solute transport on the flow field to determine the distribution of nutrients and products (e.g., Kleinstreuel & Agarwal, 1986; Salmon et al., 1988; Pillarella & Zydney, 1990). It has been shown that, under typical HFBR operating conditions, diffusion is the primary transport mechanism for low-molecular-weight nutrients and metabolites (Webster & Shuler, 1978; Piret & Cooney, 1991) while convection is the dominant mechanism for transport of macromolecular species such as growth-factor and product proteins (Taylor et al., 1994). Only recently has the effect of the osmotic pressure on the ECS hydrodynamics and protein distribution in ultrafiltration HFBRs been modelled (Patkar et al., in press; Taylor et al., 1994). A summary of all the factors which affect the convective flow in HFBRs is presented in Table 2.1. The following sections present the general assumptions used by the existing models and then more detailed descriptions with an emphasis on models based on the Krogh cylinder approximation. 2.2. General assumptions The fluid in HFBRs is assumed to be incompressible and Newtonian; body forces are neglected. Then, the fundamental equations governing the hydrodynamics become: Chapter 2: Previous modelling work 7 Table 2.1: Parameters affecting the flow in hollow-fibre bioreactors. GENERAL Design and operating parameters Fibre arrangement Membrane properties Others SPECIFIC Flow configuration (Braining, 1989; Tharakan & Chau, 1986, 1987) Flow direction Reactor orientation Geometrical design EXAMPLES AND REMARKS closed-shell (open lumen inlet and outlet) cross-flow filtration (open ECS inlet, lumen outlet) permeate suction (open both inlets and lumen outlet) dead-end filtration (open lumen inlet and ECS outlet) unidirectional flow (most cases) periodically alternated (Piret & Cooney, 1990a) horizontal (Piret & Cooney, 1990a) vertical (Patkar et al., in press) shape of lumen manifolds (Park & Chang, 1986) shape of ECS manifolds and location of ECS ports ECS and lumen inlet-outlet pressure differences (and hence the flow rates) Packing density Parallel alignment Uniformity in distribution Isotropicity Permeability, L p Nominal molecular weight cut-off Osmotic effects Cell growth conditions Membrane fouling more resistance to flow if higher (Kelsey et al., 1990) deviations from increase resistance to transverse flow and (less) to parallel flow (Kirsch & Fuchs, 1967) more resistance to flow if higher (Jackson et al., 1986) non-uniformity causes channelling (Heath et al., 1990) isotropic or anisotropic (Waterland et al., 1974) expressed as fluxXviscosity/(pressure difference) roughly, molecules with lower molecular weight will pass through the membrane (Cima, 1988) flow influenced by osmotically active species ECS porosity, hydraulic conductivity, diffusivity, viscosity, density affected by cells and proteins pore-blocking, adsorption, protein denaturation and gel-layer formation (Mulder, 1991) Chapter 2: Previous modelling work 8 continuity V - V * = 0 or V-V = 0 (2.1) momentum (Navier-Stokes) + V*-(VV*) = - - V P + ^ - V 2 V * (2.2) at P P Darcy's law V = - - V P (2.3) where V * and V are the actual and superficial velocity vectors, respectively, P is pressure, k is the Darcy permeability, p is the fluid density, /i is the fluid viscosity and t is time. Essentially, three regions are distinguished in the analysis: lumen, membrane and shell (ECS). Mass continuity applies to each of them. The Navier-Stokes equation is used for regions not treated as porous media (lumen and cell-free ECS), otherwise Darcy's law is employed (membrane and cell-packed ECS). In most models, fully developed laminar flow in the lumen is assumed and the inertial terms in Equation 2.2 are neglected because of a very small aspect ratio (fibre radius/length) and hence small radial Reynolds number (Kelsey et al., 1990; Taylor et al., 1994). This leads to the creeping-flow equation: VP = /xV2V* (2.4) For the case where proteins are present in the shell side, the solute balance equation can be written as follows: a f1 — - = V*-VC + Z)V2C + * (2.5) a t Chapter 2: Previous modelling work 9 where D is the protein difrusivity which is treated as a scalar independent of concentration. The sink/source term ¥ can include protein leakage through the hollow-fibre membrane, protein consumption or production, etc. (it is usually assumed to be zero). Owing to the reactor shape, cylindrical coordinates are used in the analysis, with the angular variation being neglected. Fibres are treated as parallel hollow cylinders distributed uniformly throughout the circular cross-section of the reactor. Only radial flow is assumed to occur in the membrane. 2.3. Major approaches to modelling In the cell models developed by Happel (1959) and Kuwabara (1959) fluid motion through an assemblage of solid cylinders was modelled by considering an equivalent system of one cylinder and a concentric fluid envelope associated with it, the ratio of solid to fluid volumes being the same as for the assemblage of cylinders. Using slightly different boundary conditions at the cell boundary (zero vorticity by Kuwabara, zero shear stress by Happel), they solved the two-dimensional creeping-motion equations for flow parallel or perpendicular to the cylinders and derived expressions for the stream function, drag and Darcy constant for the bulk flow behaviour. This simple model provides no information about the local velocity or pressure profiles in different regions of a multi-fibre reactor where the shell and lumen spaces communicate across a semipermeable membrane. Most HFBR models have been based on the assumption that the flow associated with each fibre is identical, so that a single fibre along with the fluid cylinder surrounding it is representative of the whole reactor (e.g., Kelsey et al., 1990; Pillarella & Zydney, 1990; Taylor et al., 1994). This single fibre unit is called the Krogh cylinder, in honour of Krogh (1919), who carried out early modelling work on capillaries in tissue assuming the same multi- Chapter 2: Previous modelling work 10 fibre geometry. The fibres are assumed to be arranged in a regular array with no fluid exchange between adjacent Krogh cylinders (Figure 2.1). Partial overlapping of neighbouring fibre units accounts for the void volume between them (Figure 2.1.a). Most analyses assume the closed-shell configuration, i.e. fluid enters the fibre through its lumen, then partially penetrates into the ECS in the upstream half of the fibre length and returns to the lumen in the downstream half (Figure 2.1 .b). The flow induced in the ECS by an axial pressure gradient in the lumen, in the presence of a permeable membrane, is referred to as Starling flow (Starling, 1896). Early theoretical studies on solute transport in HFBRs usually neglected the convective effects (except in the lumen) and assumed only radial diffusive transport of substrate and product in the membrane and ECS, with chemical reaction taking place in the latter (Rony, 1971; Waterland et al., 1974; Kim & Cooney, 1976; Webster & Shuler, 1978). An order-of- magnitude analysis indicates that, under typical process conditions, salts, non-electrolytes, and gases in the ECS of HFBRs are primarily transported by radial diffusion (Piret & Cooney, 1991). Kleinstreuel and Agarwal (1986) have simultaneously solved the transient convective- diffusion, Navier-Stokes, and continuity equations assuming no axial velocity in the membrane or in the "spongy matrix" region (ECS packed with biocatalyst). However, their model ignored both the dependence of hydrodynamics on mass transfer and the convective transport in the ECS. Apelblat et al. (1974) performed a theoretical analysis of a thin-walled capillary surrounded by a porous bed of tissue, with the flow in the latter described by Darcy's law. This situation is analogous to a hollow-fibre bioreactor with a densely packed cell bed in the ECS. The coupled steady-state continuity and momentum equations were solved for the lumen and surrounding tissue, with the results presented in terms of Bessel functions. Salmon et al. (1988) extended Apelblat's analysis to describe fluid flow as well as solute Chapter 2: Previous modelling work 11 a. Figure 2.1: Krogh cylinder approximation: a) fibre arrangement throughout the reactor cross-section, b) longitudinal section of a single unit with major flow paths. transport in a specially designed hollow-fibre reactor that consisted of the following four regions: lumen, inner membrane (permeable to fluid flow), cell- or enzyme-packed annulus, and outer membrane (impermeable to fluid flow) (Libicki et al, 1988). Their model included both convective and diffusive transport of either a non-reactive tracer or a solute consumed with first-order or Michaelis-Menten kinetics. The convective transport was found to have a marked effect on the reactor performance under extreme flow conditions. In some cases, the Chapter 2: Previous modelling work 12 transport of non-reacting solute could be adequately described using axially-averaged velocities, producing a simpler model whose solution required much less computational effort. Bruining (1989) presented a general description of the hydrodynamics in hollow-fibre devices. The scope of his analysis included different modes of operation (e.g., closed-shell, continuous open-shell, suction of permeate, dead-end filtration) corresponding to various applications of hollow-fibre modules. Starting from the mass and momentum balance equations, Bruining obtained expressions for the hydrostatic pressure and bypass (fraction of fluid passing through the ECS) as functions of the axial position and the dimensionless transport modulus TL (the ratio of the viscous resistance inside the fibre lumen and the permeation resistance of the membrane). Bruining's simple analysis provided no information on local velocity profiles. Another hydrodynamic model was developed by Kelsey et al. (1990) whose analysis was similar to that of Apelblat et al. (1974) except that the Navier-Stokes equation rather than Darcy's law was employed for the cell-free ECS. Lumen and shell pressures were assumed to be radially constant. Steady-state analytical solutions were obtained in terms of three dimensionless parameters: y, which describes the geometry of the hollow-fibre module, K , the membrane permeability and / , the filtration fraction (fraction of fluid leaving the device through the ECS downstream port). Pillarella and Zydney (1990) extended Kelsey's analysis to include glucose and insulin transport in a hollow-fibre bioartificial pancreas where both solutes were present at low concentrations. Because osmotic effects were unimportant in this case, the flow equations were decoupled from the substrate and product transport equations. Axial diffusive transport throughout the reactor was neglected and only radial flow in the membrane was permitted. The steady-state fluid flow profiles were evaluated analytically (Kelsey et al., 1990) while the transient convective-diffusion equations for glucose and insulin were solved numerically. The Chapter 2: Previous modelling work 13 model predictions were in good agreement with experimental data obtained by Colton et al. (1980). Taylor et al. (1994) incorporated osmotic effects and solved the two-dimensional protein transport equation coupled with a second-order ordinary differential equation for the radially- averaged lumen (or ECS) velocity. A transient solution was obtained by iterating the interdependent velocity and concentration fields at each new time step. Taylor's analysis was carried out for both single- and multi-fibre isotropic membrane HFBRs with an ECS essentially unobstructed by cells (i.e. during the start-up phase). The results confirmed the occurrence of a significant downstream polarization of ECS proteins. It was found that, at higher protein concentrations and lower recycle flow rates, the osmotic influence of the proteins could reduce the Starling flow by several orders of magnitude, thus eliminating the protein polarization problem. It was suggested that introducing a concentrated solution of inert, osmotically active, macromolecules with the inoculum would allow more rapid and uniform cell growth in HFBRs, leading to reduced start-up time and increased reactor productivity. Patkar et al. (in press) followed up Taylor's work and compared predictions of one- and two-dimensional models with experimental results, concluding that the radial variations could be neglected. Good agreement was found for experimental and theoretical transient and steady-state axial concentration profiles of bovine serum albumin (BSA) and human transferrin. At the upstream and downstream ends of the ECS some discrepancies were observed, which were believed to be due to the presence of the ECS manifolds. The influence of the flow direction switching time in the bidirectional lumen flow mode and the effect of membrane permeability on the protein distribution were also investigated. A formula for the critical protein loading necessary to ensure that the steady-state growth factor distribution would extend over the full length of the ECS was developed. Chapter 2: Previous modelling work 14 Koska (1993 a) recently investigated protein redistribution in HFBRs with a gel-packed ECS. Experimentally obtained protein concentration profiles were compared with those predicted by one- and two-dimensional models based on the Krogh cylinder approximation. The results indicated that the one-dimensional model, which required about two orders of magnitude less computational time, was sufficient to adequately duplicate the ECS protein distributions predicted by the two-dimensional model. Koska's model simulations also showed that protein polarization, a dominant feature of cell-free ECS protein transport, was reduced under cell-packed conditions, when the ECS hydraulic conductivity was lower. The governing equations of the Krogh-cylinder-based models of Kelsey et al. (1990), Taylor et al. (1994), Patkar et al. (in press), and Koska (1993a), are included in Appendix A. Chapter 3: Development of the Porous Medium Model 15 Chapter 3: DEVELOPMENT OF THE POROUS MEDIUM MODEL 3.1. Introduction Since a hollow-fibre reactor contains thousands of densely packed fibres, an attempt to describe it as a porous bed seems well justified. In the Porous Medium Model (PMM), the shell (ECS) and lumen sides are treated as interpenetrating porous regions with a continuous, spatially dependent, source/sink of fluid. This approach is analogous to the model of flow in tissue proposed by Baxter and Jain (1989). Since fluid incompressibility and Darcy's law are used to describe the hydrodynamics in both the ECS and the lumen side, the hollow-fibre membranes do not have to be distinguished as a separate region. Axial flow in the membrane is neglected, which, together with the incompressibility requirement, implies that what fluid disappears from (or appears in) the lumen must instantly appear in (or disappear from) the ECS (at the same position). The protein is assumed to be present in the ECS only, with no leakage into the fibre lumina and no build-up in the membranes. Osmotic effects are included and cause the coupling of the fluid flow and protein transport equations. Body forces (gravity) are neglected and axial symmetry of the system is assumed. In the derivation of the model equations, it is useful to introduce the concept of the representative elementary volume (REV) (Bear, 1972). The REV must be small enough to ensure continuous and smooth variations of concentration and flow properties over the length and cross-section of the HFBR. On the other hand, it must contain a sufficiently large number of fibres, so that its actual heterogeneity is not pronounced. Uniformity in fibre distribution is not essential, although it is convenient to assume that each REV contains the same number of Chapter 3: Development of the Porous Medium Model 16 R D B ^Ax. Figure 3.1: Ring-shaped representative elementary volume (REV) with thickness Ar and length Ax. fibres per unit volume. A diagram of a two-dimensional REV in cylindrical co-ordinates is shown in Figure 3.1. The PMM is the first attempt to develop a more general HFBR model, able to deal with a variety of flow configurations, with possible significant radial pressure and concentration gradients and, eventually, with non-ideal reactor design details (e.g., ECS manifolds) or the effect of gravity. In contrast to the Krogh cylinder approach, in which a multi-fibre reactor is modelled by considering a fictitious single fibre unit, the spatial domain in the PMM corresponds to the real dimensions of the HFBR cartridge. Chapter 3: Development of the Porous Medium Model 17 3.2. Hydrodynamics It is assumed here that, because all of the components of the system are incompressible, the reactor hydrodynamics are always quasi-steady. Any transient changes in the flow field are due to the time-dependent changes in protein concentrations (via osmotic pressure). Thus, at each time level, the hydrodynamics adjust instantaneously to the new concentration field. As well as being incompressible, the fluid is also assumed to be Newtonian. The steady-state continuity law applied to the ECS yields V-Vs=</> (3.1) and to the lumen, V - V L - - * (3.2) where Vs and VL are the shell and lumen superficial velocity vectors, respectively. The fluid source/sink term, <t>, is due to fluid leakage across the membrane and can be expressed as 0 = - ^ ^ ( P L - P s + n s ) (3.3) /* where Lp is the membrane hydraulic permeability, /x is the fluid viscosity, Ay is the membrane surface area per unit volume available for fluid transport, PL(x,r) and Ps(x,r) are the lumen and ECS hydrostatic pressures, respectively, and ns(x,r) is the ECS osmotic pressure. As mentioned before, it is convenient to assume that Ay is a constant independent of position and time. In the absence of any membrane fouling phenomena (Mulder, 1991), the same assumption can be made for Lp. The fluid viscosity is assumed to be independent of protein Chapter 3: Development of the Porous Medium Model 18 concentration, which, as shown by Koska (1993 a), is a reasonable approximation for the concentration range of interest here. Since the Reynolds numbers in HFBRs are usually very small (e.g., for a lumen flow of 600 cmVmin, the lumen Re=o(l); ECS Reynolds numbers can be orders of magnitude smaller), Darcy's law can be employed, thus giving a simple relationship between the local velocity components and corresponding pressure gradients. In cylindrical co-ordinates, with angular terms neglected, Darcy's law becomes V - - I f a n a rt \ d P« . , dR K K,s ~~2 + l r K.s ""a" ^ ox or (3.4) for the ECS, and ox dr (3.5) for the lumen (lx and lr are unit vectors in the axial and radial directions, respectively). The principal components of the hydraulic conductivity tensor, kx and kr, are assumed to be constant throughout each medium, although, in general, different in either direction. Moreover, since the fibres are not directly connected with one another, the lumen flow is essentially one-dimensional and krL can be set to zero, yielding 1 dP VL = - L K K , L ^ - (3-6) fl o x Section 3.4 discusses in more detail how the hydraulic conductivities are modelled. Combination of Eqs. 3.1 with 3.4 and 3.2 with 3.6, with regard to Eq. 3.3, yields the following set of coupled partial differential equations for PL(x,r) and Ps(x,r): Chapter 3: Development of the Porous Medium Model 19 - ^ ^ - = ^ ^ ( P L - P s + n s ) (3.7) ox * J t f ^ t - ^ 4 r ( p L - p . + n8) (3.8) dx For identical, straight fibres, the surface area per unit volume, Ay, can be expressed as _ 2rRLLn _ 2RLn A ~ *R*L " ~ # ~ ( 3 9 ) where n is the total number of fibres in the HFBR, R is the cartridge inner radius, L is the reactor length (i.e. the ECS length) and RL is the fibre inner radius. The actual value of Ay may be larger than that calculated from Eq. 3.9 because of fibre swelling in the liquid-filled cartridge (Patkar et al., in press). The hollow fibres are reported to assume a wavy appearance and have both their radial and axial dimensions increased by about 10%. Note that the determination of Lp should be based on the surface area Ay calculated from Eq. 3.9. Tables 5.1 and 5.3 of Chapter 5 list the numerical values of all the parameters used in the model. 3.3. Protein transport Since the ECS is assumed to be the only region that contains protein, there is just one differential equation describing the protein transport. It has the general form of Eq. 2.5, with the sink/source term SF set to zero, as protein leakage, adsorption, denaturation, production, and consumption are neglected. The equation can be derived from a protein mass balance over the representative elementary volume. Figure 3.2 shows the enlarged REV cross-section r d r 3jk 3r Chapter 3: Development of the Porous Medium Model 20 D_ N C mmmmmmmmmmmmmmmmmm Sm588^$i8SS&$$K8^8888S«888S8£S^^ sssssssssssssssssssss^^ S ^ S ^ S S S S S ^ ^ g ^ ^ sssss^^s^ss^sss^s^ SSSSS-WSSSSWSSSSSSg ̂ iSWSSSSSSMSSKSiSSS ^ S S S ^ S ^ ^ ^ ^ ^ ^ ^ fe>,^^^^^5SSSS^S^ : g^sss^sss^^^^ss^^g^^s^^ss^ss I I Figure 3.2: The REV cross-section diagram for the convective-diffusion equation. marked ABCD in Figure 3.1. The rate of protein accumulation in the REV must be balanced by the net convective and diffusive fluxes through all the boundaries, according to the following equation: 3 c e^x(r 2 2 - r i 2 )Ax— = accumulation T ^ - l f ) u c - A a x - *"(r2 2 -ii2) , 9 (u c) „ u c + A x —^—- - Dx ox (3.10) 3 c 82c -— + A X -—r- ax dx2 W face E face + 2 x r, A x ~ 9c v c - Z ) r - 2 xr2 Ax * 3 v c ^ vc + Ar , -Dr 3r 3 c d2c a r a r S face N face Here, c is the average actual concentration in the REV (the mass of protein in the part of the representative elementary volume available to the fluid), u is the superficial velocity at face W and v is the superficial velocity at face 5. Dx and Dr are the effective diffusivities of protein in Chapter 3: Development of the Porous Medium Model 21 the axial and radial directions, respectively. es is the overall porosity of the HFBR, expressed as follows: &S = ^ ECS ^ ECS (•*•*•*•) where 1 - e ^ is the fraction of the reactor volume occupied by fibres (including their porous membranes and lumina) and 1 - e ^ 5 is the fraction of the ECS occupied by cells. Thus, es represents the fraction of the reactor volume available for the ECS fluid. In the absence of cells, ss = eECS. In the present analysis, it is assumed that eECS and s'^g (and hence also es) are constants independent of time and position. A similar assumption has been made with respect to Dx and Dr. Moreover, variation of the difrusivities with protein concentration has been shown by Koska (1993a) to be insignificant for protein loadings below 100 kg/m3 and, therefore, has been neglected here. Section 3.4 will explain in more detail how the axial and radial difrusivities are modelled. It should be noted that, in writing Eq. 3.10, es has been incorporated in the expressions for Dx and Dr. Cancelling identical terms with opposite signs and dividing both sides of Eq. 3.10 by 7r(r22 - r,2) Ax yields dc _ d2c a (uc) 2 f n 3 c "l 2r2 (n d2c 9 (vc)^ es — = Dx — -—- + Dr v c + — Dr — —L \ at ox dx fi + r2 v ^ r J r i+ r2V ^ r ^ r J Since Ar is assumed to be small, r, «r2 « r , which eventually leads to e ^ 7 7 = ^~\D'1~"UC\ + -T- r A T ~ - r v c . (3.12) dt axy ox J r d r ^ dr J The above second-order, time-dependent, convective-difrusion partial differential equation is coupled with the pressure equations (3.7 and 3.8) through a relationship n s (c ) between Chapter 3: Development of the Porous Medium Model 22 osmotic pressure and concentration. In this work, bovine serum albumin (BSA) has been chosen as a model protein. Its physical and chemical properties are well described in the literature and the protein is relatively inexpensive for experimental study. The osmotic pressure of BSA can be expressed using, for instance, the following formula obtained by Vilkeretal. (1981): H,(c) = ^(j(zpcf+(2msMp)2 -2msMp + c + A2c2+A,c^ (3.13) where Rg is the gas law constant, T is the absolute temperature, Mp = 69 kg/mol is the molecular weight of BSA, Zp is the protein charge number and ms is the molar salt concentration. The virial coefficients ^ and 4? a r e functions of Zp. Refer to Table 5.1 for the numerical values of the parameters used in this relationship. The superficial velocities u and v in Eq. 3.12 are calculated from the ECS hydrostatic pressure field obtained as a solution of Eqs. 3.7 and 3.8 at each time step, i.e., 1 , 5 Ps u = Ks-r^ (3.14) H ox v = - 1 ^ ^ . (3.15) fi or The actual velocities, u* and v*, are related to u and v through es: u* = — u (3.16) v* = — v . (3.17) £* Chapter 3: Development of the Porous Medium Model 23 3.4. Modelling of the ECS and lumen hydraulic conductivities and protein effective diffusivities Hydraulic conductivities in the lumen and cell-free ECS A simple one-dimensional analysis of the laminar Krogh cylinder flow (Kelsey et al., 1990; Taylor et al., 1994) provides the following expressions for the radially-averaged actual axial velocities in the lumen and cell-free ECS: uT = - $fi dx (3.18) u, = - _1_ 8/x 4% IniRs/Rj Rl-R2 -3i£+i?j M M d^L dx (3.19) where RL is the inner fibre radius, RM is the outer fibre radius, Rg is the Krogh cylinder radius and PL and Ps are the lumen and ECS pressures (assumed radially constant). It should be pointed out that Rg is calculated based on the assumption that the sum of all Krogh cylinder volumes equals the total reactor volume, i.e., R. = -T=R \n (3.20) where n is the number of fibres and R is the cartridge radius. If the actual velocities u,* and us* are converted into superficial velocities, uL and us, Eqs. 3.18 and 3.19 will essentially be expressions of Darcy's law and will yield the Darcy permeabilities for the lumen and for the ECS. The conversion is performed in an analogous way as in Eq. 3.16 and requires knowledge of the lumen and ECS porosities: Chapter 3: Development of the Porous Medium Model 24 eL = Rl/% (3.21) e* = CECS = 1-Rl/Rl (3-22) This leads to L 8/i ^ dx v ' — _1_ 8/x us = - ^ . &\ (4Rt ln(^/^) _ , D2 ^ D 2 1 dP5 i2 I ^S - 3 i ^ + i ^ ^ a . (3.24) V P4-K * " J dx which eventually yields the following expressions for the Darcy permeabilities in the axial direction: K, - i f (3-25) Rlf . 3 _ 1 *,* = ^{-^v-^2*-^*21 (326) where ^ = 1-£ECS = RM/R-S 1S t n e fraction of the reactor volume occupied by the fibres. It is worth noting that Eq. 3.26 is identical with the prediction of Happel (1959), who used a cylindrical cell model to investigate laminar flow parallel to an array of cylinders. In the absence of cells, the flow in the extracapillary space is analogous to parallel or transverse flow through a bank of cylinders. Examples of systems having a similar geometry, which have been investigated using both experimental and theoretical methods, include tubular heat exchangers (Sangani & Acrivos, 1982a; Hwang & Yao, 1986) and fibrous porous filters (Spielman & Goren, 1968; Harrop & Stenhouse, 1969; Ethier, 1991). The dimensionless Darcy permeabilities (or hydraulic conductivities), kJR2M, which can be found in the literature Chapter 3: Development of the Porous Medium Model 25 for flow parallel or perpendicular to an assemblage of parallel cylinders, are summarized in Table 3.1. To allow comparison between different authors' predictions, some values of kxS and krS have been calculated for RM = 10_4m and <p = 0.5 and are included in the Table. It should be stressed that an analysis based on the Krogh cylinder approximation cannot provide an estimate of the ECS permeability in the radial direction. In this case, one of the expressions from Table 3.1 can be employed in the model (although, some of the values listed differ significantly and the choice of krS is rather arbitrary). Here, the simple cell-model equations due to Happel (1959) were used to calculate the ECS hydraulic conductivities in both directions. Protein diffusivities A similar difficulty arises when specifying the ECS effective protein diffusivities, Dx and Dr, which are different from the diffusivity D in a fluid unobstructed by the presence of fibres or cells. Using again the analogy to an array of parallel cylinders, Neale (1977) obtained the following relationships for these effective diffusivities: A = DEECS = D{l-<p) D = D SECS = D ^ - 2-SECS 1 + P Cell-packed ECS In this case, the hydraulic conductivities cannot be obtained from the equations listed in Table 3.1 by simply including the e'ECS factor. The reason is that the flow in the cell-packed ECS is no longer analogous to that through an assemblage of cylinders (fibres) and should (3.27) (3.28) Chapter 3: Development of the Porous Medium Model 26 Table 3.1: Theoretical expressions for the dimensionless Darcy constant k/R2M as a function of the solid fraction <p for a flow through an array of parallel cylinders (RM is the cylinder radius); k values are calculated for <p=0.5 and i ? M = 1 0 - 4 m. a. Flow parallel to the cylinders Authors Expression for kxSJR^ 2 M Value of k x.S Happel (1959) Eisenberg & Grodzinsky (1988) Kelsey et al. (1990) Taylor et al. (1994) *x,S R M J_ 4<p 3 1 , - I n a ? — + 2<p- — u> 2 2 3.4-10" l um - 1 0 ^ , 2 Spielman & Goren, (1968) vx,S R J_ 2<p V^7 Jo(VV^">" •M -r L RM J , ( V V ^ ) . J p (x) is Bessel function of order p and argument x Drummond & Tahir (1984) -^- = -^\-ln<p-\A975+2<p-^<p2 RM *<P [ 2 (equilateral triangular array) + . nxS J _ 1 7 - In<p -1.4763 +2<p--<p2 3.5-10" iUm -10 ~ , 2 4.5-10" l um -10 „ 2 -0 .051^ 4 / ( l +1.5198v?4)) (square array) Chapter 3: Development of the Porous Medium Model 27 Table 3.1 - cont. b. Flow perpendicular to the cylinders Authors Happel (1959) Kuwabara (1959) Spielman & Goren, (1968) Sangani & Acrivos, (1982a) Expression for krSJR2M Ks I Rl *<P Ks I R2M 8* Ks Rl r «2 i I " l n ^ + „»4.1 3 - I n ? — + 2co 2 [— M 4<p 2 . y/Ks' h(RMljKs) _ RM J , ( W ^ 7 ) _ Jp(x) is Bessel function of order p and argument x Ks 1 Rl 8* (equilatera 1 i -\n<p-\A9 + 2<p — ? +... il triangular array) -£f = T " (" In? " 1-476 + 2<p-1.774?2 RM &(p + 4.076 <p3 +...) (square array) Value of £,. s 2.3-10~10m2 4.8-10"10m2 • 2.0-10~10m2 7.M0 - I 0 m 2 Chapter 3: Development of the Porous Medium Model 28 Table 3.1b-cont. Authors Drummond & Tahir (1984) Eisenberg & Grodzinsky (1988) Expression for krSJR2M *r,S 1 Rl,~*<p - \n<p-1.4975 +2<p--<p2 - 0 . 7 3 9 1 / + ...) (equilateral triangular array) K « 1 + (l<p - 0J959<p2)/(l + 0.4892<p -1.6049 ?2)) (square array) Ks i R2M 8*>G [-ln*V+ij G is a function of porosity, viscosity, dielectric permittivity, surface charge, bulk fluid conductivity, double layer thickness, counter ion mobility, potential difference across double layer Value of kr s 0.6-10-10m2 4.2-10~10m2 Chapter 3: Development of the Porous Medium Model 29 rather be modelled as a flow through a bed of densely packed spheres (cells). This implies that the conductivities and diffusivities in both directions become essentially equal, as in an isotropic medium. Although the maximum packing density of rigid spheres is 74% (i.e., 26% porosity, Bear, 1972), cells are deformable and hence can be packed to higher densities. Here, it will be assumed that effective difiUsivities and Darcy permeabilities of packed beds of compressible cells can be described, by extrapolation, using relationships for a bed of spheres having the same diameter as the cells. Since neither the effective diffusivity, £)*, nor the hydraulic conductivity, k', nor the cell-packed ECS porosity, e ^ , have been determined experimentally for an HFBR packed with mammalian cells, it would be convenient if both D° and k' could be expressed as functions of s'ECS (or ss). Then, only the latter parameter (rather than all three of them) would require estimation. Accordingly, the Carman-Kozeny equation (Carman, 1937; Kozeny, 1953) can be used to express k* in terms of the cell-packed ECS porosity, i.e., *• = ( £ ^ ) 3 (3.29) where 4 / 1S the porous medium surface-area-to-volume ratio. Assuming that the cells are spheres with a 12 /xm radius, we obtain from Eq. 3.29 the values of k' =1.03-10~13 m2 and k' =4.43-10"16 m2 for ss=26% and 5%, respectively. An estimate for k' has been provided by Koska (1993a) who measured the hydraulic conductivity of an agarose gel that simulated the cell-packed environment. A range of 0.9-10-16 — 20-10~16 m2 was obtained, depending on the gel concentration and on the presence of additives (BSA). The hydraulic conductivities of natural cell aggregates of Escherichia coli packed in the ECS of hollow-fibre bioreactors were measured by Libicki et al. (1988) who obtained an approximate range of 4-10-12 — 3-10-11 m2 (depending on the Chapter 3: Development of the Porous Medium Model 30 cell volume fraction). These estimates are several orders of magnitude higher than the values calculated using either the theory developed by Sangani and Acrivos (1982b) or Eq. 3.29. The former theory yields a range of 2-10"15 — 4-10-12 m2, corresponding to e^, « 0.70 - 0.08; whereas Eq. 3.29 yields a range of 2.7-10"18 -2.6-10"13 m2, corresponding to the same s'ECS range as in Libicki et al. (1988), i.e., 0.88 - 0.075. Other relevant estimates of Darcy permeabilities are available for bacteria in filter cakes (10 ~14 - 10 ~16 m2, Humphrey et al., 1985), packed beds of red cells (7-10-15 - 3-10"18 m2, Zydney et al., 1986) as well as subcutaneous and hepatocarcinoma rat tissues (Swabb et al., 1974) or rabbit omentum tissue (Apelblat et al., 1974), both of which generally fall below 10 ~18 m2. The effective diffusivity, £>*, can be calculated from the formula derived for a bed of spheres by Neale and Nader (1973) modified by the eECS factor accounting for the presence of fibres: D' = DeECS^^. (3.30) Both Eqs. 3.29 and 3.30 were assumed to be valid for deformable spheres, i.e., for packing densities larger than the maximum for rigid spheres. 3.5. Boundary conditions In general, in order to find a solution of a given differential equation, it is necessary to specify one boundary condition per independent variable per equation order with respect to that variable. Thus, for the ECS pressure equation (Eq. 3.7), which is of second order with respect to x and of second order with respect to r, four boundary conditions will be needed. Similarly, for the lumen pressure equation (Eq. 3.8), two boundary conditions at two x- Chapter 3: Development of the Porous Medium Model 31 positions are necessary. The parabolic protein transport equation (Eq. 3.12) requires five boundary conditions: four in the spatial domain plus one initial condition in the time domain. Altogether, we shall need 11 boundary conditions for our system of equations. It should be pointed out that, in the open-shell cases, at T=R, we will have at least one extra region with different boundary conditions for x<xm and/or x>L-xm (corresponding to the upstream and downstream ECS manifolds, respectively) than for x m <x<Z-x m (see Figure 3.3). For example, instead of the no-radial-flux condition at v=R, as exists for closed-shell operation, we may have a specified inlet flux or a known pressure at the ECS upstream port, i.e. for x<x„. LUMEN INLET ECS UPSTREAM PORT n T=R ECS DOWNSTREAM PORT J l SPATIAL DOMAIN r=0 LUMEN OUTLET x=0 x=x. x-L~Xm x=L Figure 3.3: Diagram of the spatial domain boundaries in the Porous Medium Model. Table 3.2 summarizes the boundary conditions common to all flow configurations. The initial concentration field c0(x,r) can be assumed constant (e.g. after inoculation) or taken as an intermediate solution of another problem. Zero-flux conditions occur in the ECS at x=0 and at x=Z (i.e., dPs/3x = 0 and dc/dx = 0), as well as at r=R for xm<x<L-xm (i.e., 3Ps/dr = 0 and dc/dr = 0). In addition, there is symmetry of ECS pressure and concentration about the centre-line (i.e., 3Ps/3r = 0 and dc/3r = 0 at r=0). Chapter 3: Development of the Porous Medium Model 3 2 Table 3.2: Boundary conditions common to all flow configurations. # 1 2 3 4 5 6 7 8 9 X any 0 0 L L any any x„<x<L-x„ m m x„<x<L-xm I In m r any any any any any 0 0 R R Condition at t=0: c = c0(x,r) dPs/dx = 0 dc/dx = 0 3Ps/dx = 0 dc/dx = 0 3Ps/3r = 0 3c/ar = 0 dPs/3r = 0 3c/3r = 0 Meaning initial condition no axial fluid flux no axial protein flux no axial fluid flux no axial protein flux symmetry about the centre line symmetry about the centre line no radial fluid flux no radial protein flux The remaining boundary conditions depend on the combination of open or closed lumen and ECS inlets and outlets. Also, in the case of the pressure equations, it is possible to specify either a known pressure or a known flow rate (and hence velocity) at the inlet or outlet. These two options correspond to the Dirichlet (specified value) or Neumann (specified derivative) type of boundary conditions. At least one Dirichlet-type condition is needed to set the reference pressure level. In this work, the known-pressure rather than known-flow-rate Chapter 3: Development of the Porous Medium Model 33 (or pressure derivative) boundary conditions were applied whenever possible, in order to improve the convergence rate of the pressure solutions. In the case of the transport equation, a known concentration can be specified at the ECS upstream port, if open, such as during inoculation. However, the outlet concentration must be calculated from the ECS concentration field, as it cannot be specified a priori. This means a zero-derivative boundary condition must be specified at the ECS downstream port, if open, such as during harvesting. Since the ECS manifolds are fibre-free regions with virtually no resistance to the flow, it has been assumed here that, in an open-shell case, the relevant pressure or pressure derivative boundary condition at x=R is constant over the manifold length, i.e. for 0<x<x m or L-xm<x<L. The same assumption has been made with respect to concentration. For similar reasons, pressure and concentration are assumed to be constant in the ECS manifold at all tangential positions measured from the ECS port. Thus, even for open-shell operation, the axial symmetry of the reactor is preserved. Table 3.3 gives the remaining boundary conditions not listed in Table 3.2. The combination corresponding to a specific mode of operation can be found without much difficulty. For instance, the closed-shell mode corresponds to the open lumen inlet and outlet with both ECS ports closed; the inoculation phase requires that the ECS upstream port be open and the downstream one closed, while both lumen inlet and outlet are open; etc. The known inlet and outlet lumen pressures or velocities can, in general, be functions of radial position. It should be pointed out that, since the first-order Darcy's law has been employed in the Porous Medium Model, the no-slip condition at the reactor walls is not specified. The extra boundary condition(s) could be incorporated if a higher-order equation (e.g., Brinkman equation) had been used for the ECS. Chapter 3: Development of the Porous Medium Model 34 Table 3.3: Boundary conditions specific for different flow configurations. # la lb lc 2a 2b 2c 3a 3b 3c Configuration lumen inlet open lumen inlet open lumen inlet closed lumen outlet open lumen outlet open lumen outlet closed ECS upstream port open ECS upstream port open ECS upstream port closed X 0 0 0 L L L 0<x<xm m 0<x<xm fff 0 < x < x Iff r any any any any any any R R R Condition IL-Pi* apL n — k = - - c - uL0 3x kxX • ax p =p apL (i = — ii ax p =p rS L S,up ar Ks v ^ a p °=o 3r Comment known pressure known inlet velocity no axial fluid flux known pressure known outlet velocity no axial fluid flux known pressure known inlet velocity no radial fluid flux Chapter 3: Development of the Porous Medium Model 3 5 Table 3.3 - cont. # 4a 4b 5a 5b 5c 6a 6b Configuration ECS upstream port open ECS upstream port closed ECS downstream port open ECS downstream port open ECS downstream port closed ECS downstream port open ECS downstream port closed X 0 < x < * iff 0<x<xm L-x<x<L m L-x<x<L m L-x<x<L m L-x<x<L m L-xm<x<L m r R R R R R R R Condition c = c ^ = 0 dv p =p aps n a p - o dr dr ^ - 0 dr Comment known inlet concentration no radial protein flux known pressure known outlet velocity no radial fluid flux radially-constant outlet concentration no radial protein flux Chapter 4: Numerical techniques 36 Chapter 4: NUMERICAL TECHNIQUES 4.1. Solution of the pressure equations The two simultaneous second-order elliptic partial differential equations for pressures were introduced in Chapter 3. They have the following forms: r d r 3PS dr - K,s ! ^ f = LP Ar (PL - Ps + n s ) (3.7) k *.L^t= LpAv(PL-Vs+Us) (3.8) where n s (c ) is a known function of protein concentration (e.g., Eq. 3.13). In this study, the concentration field needed for the evaluation of n s (c ) is taken from the previous time step, although an iterative scheme simultaneously updating Ps, PL, and c (and hence n s (c)) at the same time level is also possible (see Section 4.3). The finite difference approach was employed to solve the above equations. To ensure fluid mass conservation, the equations were discretized by integration over a control volume, according to the scheme recommended by Patankar (1980). The control volume here is equivalent to the representative elementary volume (REV), the concept introduced in Section 3.1 of the previous chapter. A uniform grid (i.e., constant Ax and Ar) was used throughout the computations. Figure 4.1 displays the cluster of 5 adjacent cells (control volumes) that contribute to the finite difference formulation for the central cell (i,j). If the centre cell is adjacent to a boundary, the cluster consists of 4 cells, and in the case of a corner control Chapter 4: Numerical techniques 37 I r M-ij * J Ci-lj Vi-lj-l w _ 'Ui-ij+i Ui-lj 'Ui-ij-i • Vy+i Py+i. Cy+1 N Vij Cij Vy-i s Ej-i Cij-i Uij+ i" U« E Uy-l' A * Vi+lj Ci+lj Vi+lj-l U i + l j ' X Figure 4.1: Cluster of control volumes contributing to the finite difference equations for a point (ij) in the interior of the domain. Note: the lumen pressure equation has contributions from three cells only (constant j). volume, of only 3 cells. No fictitious points were explicitly used in the computations. The boundary pressure values, unless given through boundary conditions, were calculated by extrapolation from the interior of the domain assuming constant curvature of pressure profiles near the boundaries. If a zero-derivative condition occurred, a constant curvature of pressure or concentration profile near the boundary was assumed and the values of P or c at the point closest to the boundary and at its fictitious equivalent located symmetrically on the other side of the boundary were assumed equal. To avoid certain numerical difficulties (described, for instance, by Patankar, 1980), a staggered grid was used, with the velocity grid points located half-way between the pressure grid points. Thus, the radial velocity grid points were on the S and N faces, while the axial Chapter 4: Numerical techniques 38 velocity points were on the W and E faces. The program can output the velocity values either located on the faces or interpolated to the cell centres. In the derivation of the discretization equations, it was assumed that: (a) r d¥s/dr is constant over Ax (which implies a locally logarithmic radial profile of Ps); (b) dPs/dx and 3PL/dx are constant over Ar; (c) (PL — Ps +I I S ) is constant over Ar; (d) (PL — Ps) changes linearly over Ax; (e) IIS is constant over Ax. With the above assumptions, it was possible to obtain a coefficient matrix in the tridiagonal form for each spatial dimension (x and r). For example, the equation for the (i j ) point included the neighbouring unknown pressure values in the axial direction, (i-lj) and (i+1 j) , and in the radial direction, (ij-1) and (i,j+l) (see Figure 4.1). Consequently, five diagonals were filled in the coefficient matrix. The equations for PL could maintain their tridiagonal form because of the lack of a radial derivative term. However, it should be kept in mind that the lumen pressure would not, in general, be radially constant; rather, PL would vary with radial position through the source term on the right-hand side of Eq. 3.8 and possibly also because of the lumen boundary conditions. The solution of the resulting sparse system of linear equations is theoretically feasible using the standard elimination or decomposition methods. However, these direct approaches would be extremely inefficient with large numbers of grid points as well as being prone to substantial round-off accumulations. In the one-dimensional Krogh Cylinder Model (KCM) used in this (see Chapter 5) and a previous study (Taylor et al., 1994), several hundreds of grid points were often necessary to achieve the desired level of accuracy. To avoid handling enormous matrices in an attempt to obtain the solution directly, the iterative line-by-line over- relaxation method was used in this work (Anderson et al., 1984). In this technique, two of the neighbouring unknown values, either (ij-1) and (ij+1) or (i-lj) and (i+1 j) , are taken from the Chapter 4: Numerical techniques 39 previous iteration, thus making it possible to obtain a simple tridiagonal system of equations. Such a system is solved as many times as there are rows, i.e., K (j = 1,..,K), or columns, i.e., N (i = 1,..,N), which then completes one iteration cycle. The procedure is repeated until the solution no longer changes. In the case of Ps, either K equations were solved N times (sweeping over columns) or N equations were solved K times (sweeping over rows), while with PL, only the latter possibility existed. A variety of options are available in the program to perform the line-by-line over- relaxation in an optimum way. Over-relaxation parameters, a s and aL, were used to accelerate the convergence, according to the following formulae: Ps"=«sPs"+( l -«s)Ps n i r 1 (4-1) PM J=«LPM+(i-«L)pL n ;~ l (4.2) where IT is the iteration counter and l < a s , a L < 2 . Furthermore, it was possible to choose the direction of sweep (over ascending or descending row/column index), the mode of sweeping (over rows, or columns, or alternately in both directions) and to either stop iterating the pressure that has converged first or carry on the iteration loop until both Ps and PL have converged. Convergence was assumed to be attained when the following condition was satisfied: max prr _ prr-i ij rij <EPS (4.3) where EPS, the convergence criterion, was set to a level that ensured a satisfactory mass balance of the fluid (typically, 10'8 Pa for Ps and 10"6 Pa for PL). The function max() represents the maximum value of all i- and j-indexed arguments, i.e., in this case, the Chapter 4: Numerical techniques 40 maximum absolute value of the difference between the local pressures in the previous and current iterations. 4.2. Solution of the convective-diffusion equation The parabolic second-order partial differential equation for concentration, es — = — D u c + xDr r v c (3.12) dt dxy dx J r d r ^ dr J was solved using the well known Alternate Direction Implicit (ADI) method (e.g., Lapidus & Pinder, 1982; Anderson et al., 1984). As in the case of the pressure equations, integration over a representative control volume was performed to ensure protein mass conservation. The concentration grid points were located in the cell centres (Figure 4.1). The boundary values were either given explicitly through boundary conditions or extrapolated from the interior points by assuming constant curvature of the concentration profiles near the boundary (for a given row/column) with derivative boundary conditions taken into account, if necessary. In both axial and radial directions, Patankar's power-law scheme (Patankar, 1980) was employed to express the concentration at each face common to two adjacent cells in terms of the concentrations in the centres of these cells. By including the effect of the Peclet number, the scheme offers an efficient way of dealing with a wide spectrum of protein transport conditions, ranging from purely diffusive to purely convective transport. The basic formulae of this technique are presented in Appendix C. One difficulty inherent in the ADI method is that it is not conservative. The reason for this is that some of the unknowns in the discretization equations are represented by values taken from the previous time instant. According to the standard ADI procedure, the two- Chapter 4: Numerical techniques 41 dimensional parabolic equation is solved just twice at each time level (one sweep over rows and one over columns) rather than iterated as in the line-by-line over-relaxation method. The result is a numerically-originated mass imbalance that is dependent on how much the local concentration values have changed over the last time increment (At). There are two remedies for this adverse condition. First, At could be kept small enough to ensure that the local concentration changes would never exceed a desired level. Secondly, the standard ADI scheme might be extended by subsequently neglecting of the transient term (dc/dt) and solving the resulting elliptic equation using iterative line-by-line over-relaxation, which would finally yield a converged concentration solution at the given time instant. Alternatively, the line-by-line over-relaxation might be applied simultaneously to the pressure and concentration equations (see also the following section). In this work, the problem of the ADI method being non-conservative was handled by monitoring the protein mass balance and choosing a sufficiently small time step, At. The criterion for reaching the steady-state concentration distribution is analogous to Eq. 4.3, i.e., max c^ -c 1 1 - 1 •J 1J At <EPSC (4.4) where EPSC was usually set equal to 10"4 kg/(m3-s). 4.3. General computational algorithm The block diagram of the general computational algorithm employed in this study is depicted in Figure 4.2, with the pressure iteration block having the form shown in Figure 4.3. ACCF, or acceleration factor, is a factor by which the current time step size is multiplied in Chapter 4: Numerical techniques / START / 1 Initialize c0(hencells 0), p w .P U o ; t : = o PRESSURE ITERATION BLOCK 1 Calculate velocities 1 t := t +A t 1 Solve concentration equation y / s t e a d y state 7 \ . / or > 'Maximum time reached V or K N O NyMaximum number • \ o f iterations l / TYES Output ~ Calculate osmotic pressure ns(c) At:=At-ACCF m / STOP/ Figure 4.2: Block diagram of the general computational algorithm used in this study. Chapter 4: Numerical techniques 43 / START / / EXIT / Figure 4.3: The form of the pressure iteration block used in this study. Chapter 4: Numerical techniques 44 / START / NOy Solve for ff Converged^ b. / START / N O y / Solve for 1? • Solve for f| Converged P and R ? Figure 4.4: Examples of the pressure iteration block algorithm (not used in this study): a) with the lumen pressure lagged behind the ECS pressure, b) with the lumen and ECS pressures iterated until both converge. order to reduce computational times, particularly when seeking new steady-state solutions where the changes in concentration with time must obviously become very slow. Examples of other forms of the pressure iteration block are shown in Figure 4.4. The algorithm displayed in Figure 4.2 assumes that the flow change over the time scale At is so small that the pressure solutions can be lagged behind the concentration solution. Alternatively, the convective-diffiision equation could be iteratively solved together with the pressure equations (Figure 4.5). However, that results in considerably longer program execution times. Chapter 4: Numerical techniques 45 / START / 1 Initialize c0(hencells 0), Solve for PSi0,F>0; t :=0 t := t+At Calculate velocities Solve concentration equation X Calculate osmotic pressure ns(c) Solve for £ Solve for F| Converged N . ^ Q P^.Pg.andc? Steady state ? or 'Maximum time reached f \ . N O or Maximum number. . of iterations ?> YES Output / STOP / At:=At-ACCF Figure 4.5: General algorithm with mlly coupled concentration and pressure equations. Chapter 4: Numerical techniques 46 In most cases tested, the convergence was extremely slow. This was probably attributable to the fact that the ECS velocities were usually very small because of the low permeability of ultrafiltration membranes and, in the case of the cell-packed system, the low hydraulic conductivity of the ECS. This was particularly problematic in the case of closed- shell operations where the ECS hydrodynamic pressure was almost constant everywhere and only weakly linked to the known pressure on the lumen side. Thus, the optimization of certain numerical parameters in the program became essential. Prior to each run, optimum values of the over-relaxation parameters, a s and aL, were found. Any change in the program (e.g., in the number of grid points, membrane permeability, boundary conditions, etc.) unfortunately results in a change in the optimum values of a s and aL. Also, the EPS values (convergence criteria) for pressures were specified cautiously to avoid excessively strict accuracy requirements that would slow down program execution. Chapter 5: Testing and Application of the Porous Medium Model 47 Chapter 5: TESTING AND APPLICATION OF THE POROUS MEDIUM MODEL 5.1. Closed-shell operation: one-dimensional case In the case of a closed-shell operation, flows enter and leave the reactor only on the lumen side. If the upstream and downstream lumen manifold pressures are radially invariant, then it is expected that, over the multi-fibre averaging volume (REV) upon which the Porous Medium Model is based, flows in both the lumen and extracapillary spaces will be one- dimensional in the x-direction. Thus, if the radial derivative terms are neglected and Ay is calculated from Eq. 3.9, the PMM equations (Eqs. 3.7, 3.8 and 3.12) for a cell-free HFBR reduce to d2 P, 2Rr nL -Ks-^t = - ^ ( P L - P s + n s ) (5.1) d2PT 2RLnL ^ ~ d V = ~ i H L ( P L - P s + n s ) (5.2) and 9c _ d2c d(uc) , „ s dt x dx2 dx K J respectively. It is proven below that, under these circumstances and with kx L and kx s given by the same expressions as were used in the Krogh Cylinder Model, the governing equations for the one-dimensional PMM become essentially identical with those derived for the one- dimensional KCM (see also Appendix A). Chapter 5: Testing and Application of the Porous Medium Model 48 Substituting for kxL (Eq. 3.25), kxS (Eq. 3.26) and R2 = nR% (Eq. 3.20) in Eqs. 5.1 and 5.2 leads to d2Ps dx2 d2PL dx2 = = AR %R] l6RLLp SlniRjR^+ARlRl ' 2 ' ( P L - P s + n s ) -3R4S -K (V -p s +n s , ) and eventually yields d 'P , 16L l d2PL l6Lt 2 f L = ^±P dx2 Rl ( P L - P s + n s ) (5.5) where y is defined by Eq. A.7. Equations 5.4 and 5.5 are identical with Eqs. A.2 and A.1, respectively, derived by Kelsey et al. (1990) using the Krogh cylinder approximation and a cell-free ECS (s'ECS = 1), except for the osmotic pressure term included here because of the coupling of the hydrodynamics with the ECS protein concentration. One-dimensional equations for the ECS and lumen pressures having the same form as Eqs. 5.4 and 5.5 were also obtained by Koska (1993a) (Eqs. A.44 and A.34, respectively), who applied the Krogh cylinder assumption to a cell-packed ECS (i.e., £*ECS < 1). Dividing both sides of Eq. 5.3 by ss yields dc d2 c 9 (u *c) at ~ ax2 ~ ax TI = D ^ - ^ ^ (5.6) where D = Dx/es and u* is the local actual ECS velocity in the axial direction. Equation 5.6 is identical to Eq. A.32, derived by Patkar et al. (in press), or Eq. A.45, obtained by Koska Chapter 5: Testing and Application of the Porous Medium Model 49 (1993 a) (with Kc = Kd = 1 and D independent of axial position), by radially averaging the two- dimensional Krogh cylinder models for the cell-free and cell-packed cases, respectively. It should be noted that, under cell-packed conditions, the values of some of the constants in Eqs. 5.4, 5.5 and 5.6 (e.g., the ECS hydraulic conductivity and protein diffusivity) are different than in the cell-free case. To help verify the correctness of the numerical code developed for the Porous Medium Model, one-dimensional solutions generated by the PMM with 3 radial increments were compared with those obtained from the equations derived above for the closed-shell case. To that end, the ordinary differential equations (5.4 and 5.5) as well as the parabolic partial differential equation (5.6) were solved using Keller's box method (e.g., Anderson et al., 1984) with some modifications to improve the rate of convergence. The important parameters used in the two test cases described below are listed in Table 5.1. The transient changes in ECS protein concentration as a function of axial position in an Amicon HFBR were monitored until a steady-state polarization was achieved. At each time instant, the local concentration values obtained by both one-dimensional models were identical within the desired accuracy (10"2 kg/m3). The time at which steady state is attained depends on the convergence criterion which was less stringent in the PMM case (because of longer program execution times). The initial concentration, Co, was assumed uniform over the length of the reactor. After approximately two hours of real time operation, the axial concentration profile remained unchanged and was indistinguishable from that representing steady state. The family of curves shown in Figure 5.1 illustrates the progression of downstream polarization of protein under predominantly convective transport conditions (Piret & Cooney, 1990a). At steady state, the osmotic pressure exerted by the bulk of the protein accumulated downstream is so high that the total ECS pressure (Ps - n s ) counteracts the pressure on the lumen side. Consequently, no fluid is exchanged between the ECS and fibre lumina and, as a result, the flow is completely shut down in this part of the ECS. Chapter 5: Testing and Application of the Porous Medium Model 5 0 Table 5.1: Parameters used in the comparative study of the PMM and KCM (Amicon HFBR). i?L = l . l -10" 4 m RM = 1.9 • 10 ~4 m i?s = 2.7-10" 4m Z = 0.2m £,=1.25-10 "13m n = 5000 jt = 6.95-10"4kg/m/s Co = 10 kg/m3 APL = 4572.2 Pa (corresponding to the radially- T=310K ms = 150mol/m3 Z, = -20.4 Mp = 69 kg/mol ^ = -5.625-10~ 4-2.4M0~ 4 4 = 2.95-10"s-1.051-10~6Z;1 D = 10"10 m2/s JD=10-4-10-10m2/s ^ i = 2.510-10"10m2 Ai,s=1.28610-9m2 lumen radius outer fibre radius Krogh cylinder radius HFBR length membrane permeability number of fibres viscosity of water at 37°C (310 K) initial concentration lumen pressure drop over the length L averaged lumen inlet velocity ML*0 = 0.05 m/s) Zp -3.664-10 +1.762-10 "7 absolute temperature parameters in the relationship between osmotic pressure and concentration (Eq.3.13) " 5 Z] = - 0.0108942 Zj = 0.0001243 difrusivity, case 1 diffusivity, case 2 axial permeabilities in the lumen and ECS (from Eqs. 3.25 and 3.26) Chapter 5: Testing and Application of the Porous Medium Model 51 Figure 5.1: Radially-averaged protein concentration in the ECS as a function of axial position and time. Another comparative study between the one-dimensional PMM and KCM was carried out for different hypothetical diffusivities (D) of the ECS protein. D was varied from 10"4 m2/s to 10'10 m2/s (although protein diffusivities greater than 5-10"10 m2/s are not realistic), corresponding to a transition from a diffusion-dominated regime to a convection-dominated regime. Both models produced exactly the same steady-state concentration curves, shown in Figure 5.2. The observed differences in the maximum slope of these profiles are due to the different relative strengths of convection and diffusion in each case. For instance, at D = 10"10 m2/s, a very steep concentration gradient (at x « 16 cm) must develop at steady state in order Chapter 5: Testing and Application of the Porous Medium Model 52 that the diffusive transfer can locally balance the convective transport of the ECS protein. Table 5.2 summarizes the convergence properties of the two models for this test case. As can be seen from the table, the time needed to reach steady state increased dramatically when the protein diffusivity was decreased and convective transport became dominant. Further decreases in D (below 10"6 m2/s) produced no significant changes in the time needed to achieve steady state. It should be noted that, in the case of very low diffiisivities, it was necessary to use more axial grid points than usually because of the steep concentration gradients obtained under convection-dominant transport of protein. 0 2 4 6 8 10 12 14 16 18 20 Axial position (err) Figure 5.2: Steady-state protein concentration profiles for different hypothetical diffusivities. Chapter 5: Testing and Application of the Porous Medium Model 5 3 Table 5.2: Summary of the steady-state runs with different hypothetical diffusivities for the PMM and KCM (Amicon HFBR). D (m2/s) IO-4 io-5 10^ IO"7 IO-8 IO"10 Number of grid points PMM 9x500 9x500 9x500 9x500 3x900 3x900 KCM 500 500 500 500 900 900 Initial time step (s) PMM & KCM 1.0 1.0 1.0 1.0 1.0 1.0 Time acceleration factor, ACCF PMM & KCM 1.001672 1 1.001672 1.001672 1.001672 1.001672 Time to steady state PMM 7 min 17 s 1 h 37 min 5 h 07 min 4 h 02 min 3 h 39 min 3 h 23 min KCM 5 min 36 s 1 h 07 min 3 h 08 min 2 h 13 min 2 h 29 min 2 h 57 min The aim of the study outlined in this section has been to demonstrate that the Porous Medium Model will give correct predictions for one-dimensional closed-shell operations, both in the convective and diffusive regimes. Perhaps of greater significance is the fact that, in the absence of the macroscopic radial gradients, the PMM actually reduces to the one-dimensional Krogh Cylinder Model, which has already been shown (Patkar et al., in press) to yield excellent agreement with experimental results for both transient and steady-state closed-shell HFBR operations. The following sections of Chapter 5 describe further studies aimed at testing the model in two-dimensional situations imposed either by allowing a radial pressure variation in the lumen manifolds (Section 5.2) or by opening one or both of the ECS ports (Sections 5.3 - 5.6). Chapter 5: Testing and Application of the Porous Medium Model 54 5.2. Closed-shell operation: two-dimensional case with inlet and outlet radial lumen pressure gradients In their hydrodynamic study of hollow-fibre devices operated in the closed-shell mode, Park and Chang (1986) found that, under some circumstances, significant negative and positive radial gradients of hydrostatic pressure could develop in the upstream and downstream lumen manifolds, respectively. These lumen manifold pressure variations led to non-uniformities in flow through the fibre lumina. Although laboratory measurements by Koska (1993 a) and Patkar et al. (in press) demonstrated that these effects were negligible for the Amicon and Gambro modules they investigated, artificially imposed radial lumen pressure gradients at the inlet and outlet of an HFBR can produce two-dimensional variations of the fluid flow and protein concentration fields in the ECS and thus provide an interesting case study to further test the Porous Medium Model. Accordingly, Dirichlet-type boundary conditions for the lumen pressure at x = 0 and x = L were specified as follows: PL.O (r) = APL(1- 0.5 r/R) (5.7a) PL,N(r) = APL-0.5vlR (5.7b) where APL = 4572.2 Pa is the axial lumen pressure drop along the centre-line of the HFBR cartridge (and is the same value as was used in the one-dimensional case discussed in Section 5.1). Equations 5.7 imply that there is no axial lumen pressure drop at the cartridge wall (r = R) and that the lumen pressure drop over the cartridge radius at either manifold equals half of that along the centre-line over the length/,, i.e., 2286.1 Pa. This strong radial variation of both PL,O and PL,N is not likely to occur in any real operation involving hollow-fibre devices, but it was purposefully imposed in order to magnify the resulting two-dimensional effects. The parameters used in this study are once again those listed in Table 5.1 (Amicon HFBR), Chapter 5: Testing and Application of the Porous Medium Model 55 except for th&Lp value which was increased from 1.2510"13 m to 510"13 m in order to enhance the ECS flow and thus improve the rate of numerical convergence. The cartridge radius, R, corresponding to the Krogh cylinder radius given in Table 5.1 equals 0.0191 m (see Eq. 3.20). The ECS hydraulic conductivity in the radial direction, kr,s, was calculated from Happel's formula (see Table 3.1b) and had the value of 0.8808-10"9 m2, which is of the same order as the kAs value in Table 5.1. A protein difiusivity of 10"10 m2/s was used and the initially uniform ECS protein concentrations, c0, of 10 and 20 kg/m3 were chosen. The main purpose of the study was to investigate the steady-state ECS and lumen hydrodynamics and the ECS protein concentration field in the presence of strong radial gradients. In the absence of protein, as expected, the lumen and ECS flow fields display a fore-and- aft symmetry (the two-dimensional distributions of the lumen and ECS velocity components are symmetric about the half-length of the reactor). The radial velocities in the ECS are positive in the upstream half and negative in the downstream half of the HFBR, i.e., the ECS fluid travels radially outward in the former region and towards the centre-line in the latter. The ratio of magnitudes of the local axial-to-radial ECS velocity components is close to 10, which, approximately, equals the reciprocal aspect ratio of the HFBR, LIR. An important consequence of the imposition of radial pressure gradients is a decrease in the average magnitude of the ECS axial velocity (by a factor of 3, compared to the situation where no radial gradients are present) and, hence, in the magnitude of the ECS convective flow. Figure 5.3 shows how the steady-state distribution of the ECS protein is affected by radial pressure gradients in the lumen manifolds (all the concentration distributions shown here were obtained using the same two-dimensional code; Figures 5.3a and 5.3c were produced with APL = 4572.2 Pa, but with no lumen manifold radial pressure gradients). In the transient phase, protein accumulation generally follows the direction of the ECS flow until, at steady state, the total ECS pressures (Ps - n s ) locally balance the lumen pressures and the Chapter 5: Testing and Application of the Porous Medium Model 56 a. S ' T '0N 0 R ^ S C. 0 S / ^ V ^ 7 F ^ 5 Figure 5.3: Steady-state ECS protein concentration in the absence (a,c) and presence (b,d) of radial pressure gradients; (a,b) c0 = 10 kg/m3, (c,d) c0 = 20 kg/m3. Chapter 5: Testing and Application of the Porous Medium Model 5 7 protein distribution becomes a consequence of the distribution of pressure on the lumen side. If radial pressure gradients are imposed, the protein is distributed over a larger portion of the ECS volume than when only axial gradients exist, reflecting the much smaller range of lumen pressures that must be countered in the peripheral region of the reactor. Thus, in this case, the flow is practically shut down over a larger portion of the extracapillary space. In each of the four cases displayed in Figure 5.3, once steady state is reached, the ECS fluid flow is essentially restricted to the protein-free upstream region. The maxima of the axial and radial velocity components are located, approximately, at half of the maximum axial and radial positions of this region's boundaries, respectively. The smaller the region, the less fluid passes through the ECS. The time needed to reach steady state varied from about 5 h (Figure 5.3c) to about 7 h (Figure 5.3b). In general, the transient phase was somewhat longer if radial gradients were present or if the total amount of ECS protein was lower. Increasing the protein loading from 10 to 20 kg/m3 also increased the rate of numerical convergence to steady state. 5.3. Membrane permeability determination The permeability Lp of hollow-fibre membranes is determined by measuring the transmembrane fluid flux at a known pressure drop across the membrane. Typically, the fluid of known viscosity /x enters the HFBR through the upstream lumen manifold and exits through the downstream ECS port (Figure 5.4). The lumen outlet and the upstream ECS port are closed. The pressure drop recorded is actually the difference between the pressures at the lumen inlet and the ECS outlet, i.e., AP = P -P LXr rL,0 rS, (5.8) Chapter 5: Testing and Application of the Porous Medium Model 58 To determine Lp, it is usually assumed that AP is the pressure drop across the membrane only. This assumption is reasonable for sufficiently low membrane permeabilities, as the membrane imposes the dominant resistance to the flow. With higher Lp values, the contribution of the membrane resistance to the total pressure drop (Eq. 5.8) is lower, so that the assumption will begin to break down. Figure 5.4: Flow configuration for Lp determination. By definition, L - liQ P AAP„ (5.9) where Q is the volumetric flow through the hollow-fibre device, APm is the pressure drop across the membrane and A is the total surface area of the membrane. If one assumes APm = AP, then the apparent permeability, Lp,app = ft Q/(A AP) « Lp, can be found from the slope of the linear relationship between Q and AP using Eq. 5.9. Table 5.3 includes the important parameters used in this study, while Table 5.4 lists the Q and AP values determined in the laboratory for a Gambro HFBR (Koska, 1993b) (the Gambro HFBR was henceforth used in all the studies described in this thesis). The flow rates were found by measuring the volume of water passing through the device over a known time, while AP was determined as Chapter 5: Testing and Application of the Porous Medium Model 59 the difference between the lumen inlet pressure head and the ECS outlet pressure measured using a water manometer. The membrane surface area, A, has been estimated as equal to 1.5 m2, corresponding to an approximately 10% increase in fibre dimensions due to their swelling. If it is assumed n = 0.001139 Pas, a linear fit of Eq. 5.9 to the data in Table 5.4 yields Lp>app = 6.1810'15 m. (The permeability value of 6.40-10"15 m, as originally determined by Koska (1993b) using the same data, had been calculated incorrectly.) The flow rates also can be calculated using the Porous Medium Model once the input parameters, Lp, AP, p and A, have been specified. The Q = j{AP) relationship obtained numerically withZp = 6.18-10"15 m, fi = 0.001139 Pas, A = 1.5 m2, and AP ranging from 0 to 12555 Pa, is plotted as the solid line in Figure 5.5. The observed good agreement of model predictions with experimental data is attributable to the fact that the membrane permeability here is sufficiently low and, consequently, the assumption APm « AP is reasonable. Table 5.3: Important parameters used in the membrane permeability determination study. HFBR length, L 0.215 m Lumen radius, RL 1.15 • 10"4 m Outer fibre radius, RM 1.25-10"4 m Krogh cylinder radius, Rs 1.75 • 10"4 m ECS manifold axial length, xm 0.024 m Number of fibres, n 8123 Viscosity, fi 0.001139Pas (water, 15°C) ECS axial permeability, k^s 4.773-10'10 m2 (Eq. 3.26) ECS radial permeability, KiS 3.247-10"10 m2 (Happel, see Table 3.1) Lumen axial permeability, kXiL 7.15910"10 m2 (Eq. 3.25) Chapter 5: Testing and Application of the Porous Medium Model 60 For higher membrane permeabilities, a larger fraction of the total pressure drop between the reactor inlet and outlet occurs within the fibre lumina and ECS and, hence, local values of APm will fall below AP. Thus, the assumption of APm » AP in Eq. 5.9 yields a value of the apparent membrane permeability, Lp,app, which is less than the true value, Lp. This effect is CO H CO O « o i M— o fa <D E 10- 9 - 8 - 7 - 6 - 5 - 4 - 3 - 2 - 1 - 0 H i I i i ! • EXPERIMENT PMM, Lp= 6.18*10"15 m r /.. 1 / - J 1 | j J • j 1 - - j . 1 ' • \ 1f - 9 - 8 - 7 - 6 - 5 - 4 - 3 - 2 - 1 - 0 2000 4000 6000 8000 10000 12000 Pressure drop (Pa) Figure 5.5: Volumetric flow as a function of AP: comparison of the Porous Medium Model with experiment. Chapter 5: Testing and Application of the Porous Medium Model 61 Table 5.4: Experimentally determined volumetric flows, Q, and pressure drops, AP, for a Gambro HFBR (Koska, 1993b). AP, Pa 12555 9470 9000 8226 5445 Q • 108, m3/s 10.33 7.75 7.25 6.58 4.36 observed in Figure 5.6, where Lp>app is plotted versus Lp, and in Figure 5.7, where the ratio LpiapplLp is plotted as a function of Lp. The two figures were created by first using the PMM to determine the flow rate for a given Lp and AP and then substituting this value of Q, along with AP, into Eq. 5.9 to obtain Lp_app. For Lp less than about 10"13 m, the membrane imposes most of the resistance to the flow and LPtaPp « Lp. For Lp greater than about 10"10 m, the contribution of the membrane to the total resistance to the flow is negligible and Lp>app no longer depends on Lp. The results presented in Figures 5.6 and 5.7 were obtained for AP = 10000 Pa but, because velocities and pressure gradients are linearly related in Darcy's law, they are independent of the pressure drop. Graphs such as that shown in Figure 5.7 can be used as correction plots to obtain a better estimate of the membrane permeability determined by measuring the fluid fluxes and the corresponding pressure drops. For instance, if the Lp value of the Gambro hollow-fibre membrane were 1.25-10"13 m (i.e., the value measured for the Amicon HFBR, Table 5.1), Eq. 5.9 used with the assumption APm =AP would yield a permeability value that is approximately 5% too low. This 5% difference would probably be observable as a Chapter 5: Testing and Application of the Porous Medium Model 62 discrepancy between the experimental and theoretical predictions of the flow rate as a function of AP (Figure 5.5). It should be noted that the relationship between Lp%app and Lp, predicted by the PMM, depends on the geometry of the hollow-fibre module. Therefore, the plots displayed in Figures 5.6 and 5.7 are essentially valid only for the Gambro FJFBR investigated here, although analogous graphs could easily be obtained for any other hollow-fibre systems. The present version of the PMM does not account for pressure losses in the lumen and ECS manifolds and in the associated tubing connections. Thus, unless the pressure taps are 1E-10 1 Q- 1E-11 n X» CD E 1E-12 8. a> e E •£ 1E-13 to Q. s 1E-14 1E-14 1E-13 1E-12 1E-11 1E-10 Actual membrane permeability, L (m) 1E-9 Figure 5.6: The apparent (calculated) membrane permeability, Lp,app, versus the actual (input) membrane permeability, Lp: prediction of the Porous Medium Model. Chapter 5: Testing and Application of the Porous Medium Model 63 connected directly to the inlet lumen and outlet ECS manifolds, even the Lp value predicted by the PMM may not be quite correct, although it should be more reliable than that calculated fromEq. 5.9. ! CD Q. <D C 2 . a E o> E O c ca a. a. ca "5 o CO a: 0 1 - n 01 - 0.001 - l 11 1E-14 1E-13 1E-12 1E-11 1E-10 Actual membrane permeability, Lp (m) 1E-9 Figure 5.7: The ratio of LPiapplLp versus the actual membrane permeability, Lp: prediction of the Porous Medium Model. Chapter 5: Testing and Application of the Porous Medium Model 64 5.4. Filtration hydrodynamics: comparison of the Porous Medium Model with the Krogh Cylinder Model A flow diagram for hollow-fibre filtration is shown in Figure 5.8. Figure 5.8: Flow diagram for the partial filtration mode {QUout >0) and full filtration mode (&.„,=<>). The Krogh Cylinder Model formulation presented by Kelsey et al. (1990) (see Appendix A) accounts not only for closed-shell operation but also for the filtration mode (where the downstream ECS port is open), by introducing the filtration fraction,/ defined by Eq. A.8 (Kelsey et al. referred to /as the ultrafiltration fraction). In the full filtration mode,/= 1; in partial filtration, 0 < / < 1; and in the closed-shell case,/ = 0. One potentially significant feature of the Kelsey filtration model is that it admits an ECS outflow only parallel to the fibre through the concentric annulus of the Krogh cylinder. In reality, the ECS outflow leaves radially through a circumferential manifold at the periphery of the fibre bundle. For the sake of comparison, Kelsey's equations for the lumen and ECS pressures as functions of the axial position, PL(x)andPs(x), as well as the equations for QUin,QUout andQSw0Ut, have been re- written in terms ofPL0,PLN andPs>1, rather than/ The modified equations are presented in Appendix B. Chapter 5: Testing and Application of the Porous Medium Model 65 The study of the Lpapp dependence on Lp, as defined and described in Section 5.3, has been repeated here with the aim of comparing the predictions of the Krogh Cylinder Model with those of the Porous Medium Model. The latter describes more realistically the macroscopic radial flows created by the circumferential ECS manifold. The resulting two curves (not plotted here) show qualitative and quantitative similarity to each other and are virtually indistinguishable forLp < 10"11 m. The accurate asymptotic value oiLpapp atLp-* oo is difficult to predict in the PMM, but can readily be evaluated for the KCM curve if the modified Kelsey equations presented in Appendix B are used. For/= 1, we find nirR4 lim a ,, = —r-fil + yWs,* -PL.0) (5.10) which, for Ps dn —PL0 =10000 Pa and withZ, RL, n and y values corresponding to the Gambro HFBR (see Table 5.3 and Appendix B), yields the critical (and maximum) flow rate of 3.82247-10"5 m3/s. Inserting this value into Eq. 5.9 along withal = 1.2619 m2, corresponding to the dry dimensions of the hollow fibres, yields the asymptotic value of Lp = 3.45-10'12 m. Figure 5.9 compares the inlet volumetric flow rates predicted by both models for a wide range of Lp values in the partial and full filtration modes. The curves were obtained with the pressure drops PSjn - PL,o = 10000 Pa for/= 1 and Ps,d„ - PL,o = PL,N - PL,O = AP= 10000 Pa for 0 < / < 1 (it should be noted that, in the latter case, with the pressure drops fixed in this way,/increases as Lp is increased). The ECS and lumen hydraulic conductivity values used in the PMM are listed in Table 5.3. The Krogh Cylinder Model predicts that, at high membrane permeabilities, the inlet flow rates in both filtration modes should approach the same asymptotic value of 3.82247-10"5 m3/s, which can easily be obtained for/= 1 (as calculated above using Eq. 5.10) as well as for 0 < /< 1. In the latter case, Hm QL,in - - ^ [ K / U -PL,O)+(PL.» -pL,0)h-irf-a+T)^, (511) Chapter 5: Testing and Application of the Porous Medium Model 66 which is identical to Eq. 5.10. The convergence of the two KCM curves in Figure 5.9 can also be intuitively deduced by considering the fact that the transmembrane pressure drop is close to zero at high Lp values and that, since no radial flow in the shell and lumen sides is allowed, the fluid has to travel over the same distance L, independent of whether it flows in the ECS or in the fibre lumina. Since the inlet-outlet pressure differences are the same for both/= 1 and 0 < / < 1 and since Ps,dn = PL,N in the latter case, then the limit of the inlet flow rate at high membrane permeabilities must be identical in both filtration modes. 4.5- 4.0- 3.5 t 3.0- « C C? 2.5- ST £ 2.0- -g 1 1.5- 1 ,0. 0.5 0.0 0 < f < 1 I I I I l l l l - I—I I M i l l l l l l l l l l | KROGH CYLINDER MODEL POROUS MEDIUM MODEL Tit] 1—I I I I I I 1E-14 1E-13 1E-12 1E-11 1E-10 1E-9 Membrane permeability (m) 1E-8 4.5 2.5 2.0 1.5 1.0 •0.5 0.0 1E-7 Figure 5.9: The lumen inlet volumetric flow as a function of membrane permeability in the partial and full filtration modes: comparison of the two models. Chapter 5: Testing and Application of the Porous Medium Model 67 Although it may seem that the two PMM curves in Figure 5.9 overlap at high membrane permeabilities, they do not, in fact, approach exactly the same limiting flow rate because k^s, krS and k^i have different values and also because the distances the fluid must travel before exiting the HFBR are different in the two cases. If/< 1, then some of the fluid leaves axially through the fibre lumina at x = L, while if /= 1, then all of the fluid leaves radially over the length of the ECS manifold. Numerical simulations have shown that the value selected for kr,s has only a weak effect on Qi,irh particularly in the partial filtration mode. Since the values of kx,s, Ks and kKL used here are of the same order of magnitude, no significant difference is noticeable between the two PMM inlet flow rate curves at high membrane permeabilities. The limit of the inlet flow rate at low membrane permeabilities for 0 < / < 1 in the Krogh Cylinder Model can be calculated as Km GU = " T 7 ( ^ -^.o) (5-12) yielding the value of 2.27827-10"5 m3/s. This residual flow rate is identical to that in the closed-shell mode and corresponds to no fluid passing from the lumen side to the ECS of the HFBR. The same result can be obtained from the standard Hagen-Poiseuille equation for n parallel cylinders of radius RL (White, 1991). The outlet ECS and lumen flow rates as functions of the membrane permeability in the partial filtration mode are plotted in Figure 5.10. The QL,out curve in the KCM approaches the same asymptotic value of 2.27827-10"5 m3/s at both low and high Lp values and displays a remarkable minimum at Lp « 10"11 m. If the membrane permeability is very low, practically all the fluid travels downstream inside the hollow fibres and QL,OUI equals the Hagen-Poiseuille value. If Lp is very high, the presence of the membrane does not affect the flow and Quout can again be obtained from the Hagen-Poseuille equation. In the intermediate region, as the Lp Chapter 5: Testing and Application of the Porous Medium Model 68 1E-14 1E-13 1E-12 1E-11 1E-10 1E-9 Membrane permeability (m) 1E-8 1E-7 Figure 5.10: The ECS and lumen outlet flow rates, Qs,out and Quout, as functions of membrane permeability in the partial filtration mode: comparison of the two models. value increases, an increasing amount of fluid passes across the membrane to the ECS, which causes a temporary decrease in the lumen flow until, eventually (at Lp > 10"8 m), the hydraulic throughput capacities of both the lumen and the shell sides reach their saturation points. The shapes of the PMM curves as well as the limiting values of QUin,QUaut and Qs OTf as Lp -» oo in this case cannot be found analytically, as the flow rates here are unknown functions of not only Lp but also the ECS permeabilities and the surface areas of the outlet regions. According to Darcy's law, the rate of flow through a porous medium depends on the hydraulic conductivity of the medium, on the pressure drop per unit length and on fluid Chapter 5: Testing and Application of the Porous Medium Model 69 viscosity. Because of the finite axial dimension of the ECS manifolds in the Porous Medium Model, a portion of the fluid travels in the hollow-fibre device over a distance shorter than L. Since the outlet ECS pressure, Ps,dn, is assumed constant over the axial length of the manifold, this translates into higher inlet flow rates obtained with the PMM than those predicted by the Krogh Cylinder Model (Figure 5.9). The considerable differences, at high Lp values, in the lumen and ECS outlet flow curves pertaining to the two different models (Figure 5.10) result mostly from the fact that no radial ECS flow is allowed in the KCM and that the outflow surface area in the PMM is about 6.2 times as large as that in the other model. Therefore, the rate of fluid discharge from the ECS becomes so high that only 17% of the inlet flow exits from the fibre lumina and Qi,in never reaches the Hagen-Poiseuille value at high membrane permeabilities, as was the case in the Krogh Cylinder Model. The PMM profiles in Figure 5.10 are slightly distorted in the range 10"12 m < Lp < 10"10 m, which can probably be ascribed to the same effect that is responsible for the minimum in the KCM outlet lumen flow curve (discussed above). Although an axial outflow from the ECS in the filtration mode is a fictitious concept, it has been enforced by the one-dimensional restriction of the Krogh Cylinder Model. All of the investigations described above and presented in Figures 5.9 and 5.10 were repeated using the PMM equations with boundary conditions changed so as to allow for axial rather than radial outflow from the ECS. The results were then identical for both models, since, in this case, there existed no macroscopic radial gradients and the PMM reduced to the KCM. As can be seen from the plots presented here, the differences between the PMM and the KCM become noticeable only when Lp exceeds, approximately, 10'12 m. Thus, the Krogh Cylinder Model will yield acceptable hydrodynamic predictions for most open-shell situations of practical interest, primarily because the ECS outlet flow is controlled entirely by the membrane resistance. Chapter 5: Testing and Application of the Porous Medium Model 70 5.5. Inoculation and relaxation In the inoculation phase of HFBR operation, the cell inoculum in a solution containing high-molecular-weight growth factors, is introduced to the ECS through its upstream port with the displaced fluid leaving by the outlet lumen port. The flow configuration for this case is shown in Figure 5.11. fi» -K- i W*S.up # HFBR *^L,out Pu* Figure 5.11: Flow diagram for inoculation phase of HFBR operation. The inoculum was assumed to be in a solution of bovine serum albumin (BSA) and the mammalian cell concentrations sufficiently low that their influence on the fluid flow and protein transport could be neglected. Thus, the ECS was assumed to be cell-free. The aim of this study was to trace the ECS protein redistribution with time for the following two cases: (1) cin = 5 kg/m3 , 60 min total inoculation time, PSup -PUN = 3000 Pa; (2) cin = 50 kg/m3 , 6 min total inoculation time, PSup -PUN = 3000 Pa, where cin is the protein concentration in the inoculum solution. Table 5.5 summarizes the parameters used in the inoculation study. The total inoculation time in the second case was ten-fold lower in order that the final protein content of the ECS be approximately the same as in the first case. These amounts were not exactly the same in both cases because the inlet and outlet boundary conditions were imposed through the fixed pressures, PSup andPAW, rather Chapter 5: Testing and Application of the Porous Medium Model 71 than through fixed flow rates, Q^in = QUout. Because of the presence of osmotically-active proteins, these flow rates, although always equal, decreased slightly with time. With the membrane permeability of 6.4-10'15 m, the above pressure differences result in a flow rate of approximately 1.5 cm3/min. Table 5.6 compares the two inoculation tests emphasizing the effect of the osmotic pressure. Variation of the flow rate with time in cases #1 and #3 is almost exactly linear. Protein concentration in the inoculum, cin Pressure head difference, P&up -PUN Membrane permeability, Lp HFBR length, L HFBR radius, R Lumen radius, RL Outer fibre radius, RM ECS manifold axial length, xm Number of fibres, n Diffusivity, D Temperature, T Viscosity, fi ECS axial permeability, k^s ECS radial permeability, kriS Lumen axial permeability, k^i ( l ) 5 k g / m \ (2)50kg/m3 3000 Pa 6.4-10-15m 0.215 m 0.01575 m 1.15-10"*m l ^ - l O ^ m 0.024 m 8123 10'10m2/s 288K(15°C) 0.001139 Pa-s 4.773- 10-10m2 3.247-10"10m2 7.159-10-10m2 (water, 15°C) (Eq. 3.26) (Happel, see Table 3.1) (Eq. 3.25) Table 5.5: Summary of parameters used in the inoculation study (Gambro HFBR). Chapter 5: Testing and Application of the Porous Medium Model 72 Table 5.6: The effect of osmotic pressure in the inoculation tests. # 1 2 3 4 Cin (kg/m3) 5 5 50 50 Duration 60min 60min 6min 6 min Osmotic effects included yes no yes no Total fluid passed through HFBR (cm3) 88.04 90.79 8.76 9.08 Average final concentration, CAVG (kg/m3) 5.38 5.55 5.35 5.55 Initial flow rate, £M08 (m3/s) 2.52 2.52 2.52 2.52 Final flow rate, £M08 (m3/s) 2.37 2.52 2.34 2.52 In both cases the inoculation phase was followed by a 20-hour-long relaxation in which all the inlet and outlet ports were closed. In the relaxation phase, the ECS protein continues to redistribute owing to local concentration and osmotic pressure gradients. The start-up of an HFBR would normally not include a relaxation phase of more than 1 h. The 20 h period was used here to explore the time scales of the inoculum protein redistribution in the absence of lumen recycle flow. Since uniform distribution of the growth-factor proteins over the volume of extracapillary space is important to the subsequent cell growth phase of reactor operation, the main focus of the study presented in this section was to look at the degree of spatial uniformity of the ECS protein concentration at the end of the inoculation and relaxation phases. Figure 5.12 shows how the ECS concentration contours vary with time over 1 h of inoculation in the case of the low inlet concentration, cin = 5 kg/m3. The local radial Peclet Chapter 5: Testing and Application of the Porous Medium Model 73 number at the inlet is of the order of 50 - 80 (with the average inlet radial velocity of the order of 10'5 m/s and 18 uniformly spaced radial grid points), which indicates dominance of the convective transport of protein. The observed dispersion of the concentration front is believed to be of numerical origin and it has been found that the extent of this "wash-out" zone decreases when more grid points are used. Because some of the fluid passes into the lumina leaving the protein behind, the front moves downstream with axially and radially decreasing speed, while more and more protein is being carried towards the front from the inlet zone. This results in two visible regions of maximal and increasing concentration; for example, after 60 min of inoculation, one region is near the upstream end at r « R/3 while the other extends from r » RI2 up to the cartridge wall (r = R) at x « 2/3 L. The trends observed in Figure 5.12 become even more extreme when the inoculation is continued for a few additional hours (not shown here). The maximum local concentrations monotonically increase while the protein front velocity decreases and approaches zero at the downstream boundary. This results in a highly non-uniform concentration distribution; even though the average ECS protein concentration is 5 kg/m3, the local concentrations can reach 50 kg/m3 and more, while a significant fraction of the ECS volume remains free of protein. Figure 5.13 displays the development of the concentration contours during inoculation with a 50 kg/m3 protein solution. After 6 min, the ECS contains approximately the same amount of protein as after 60 min of inoculation with a 5 kg/m3 solution. In this case, the extent of ECS penetration by the inoculum is much smaller while the local concentrations and osmotic pressures are much higher. The temporally increasing local osmotic pressures in the ECS cause a decrease in the transmembrane pressure difference and hence a decrease in the flow rate. Consequently, after 6 min of inoculation with c,„ = 50 kg/m3, the average ECS protein concentration is slightly less than after 1 h of inoculation with cin = 5 kg/m3 (Table 5.6). Chapter 5: Testing and Application of the Porous Medium Model 74 z o o a. < Q < DC n -*J \^- ECS INLET FLOW i • / / • / / / / / • ' I i ! / \ ////! I ! ! . . : i I / 1111 I • ' \ \ ' ^ 6 - ^ / ' / / / / ' x̂vS. ; ^/^p/^ , V ^ O 1 — 3 — ^ - \ > + 2-*rT 1 " ' " ^ y / / / / \ *y\S^jy / i yr\yy^y • . . - . . • . \ ^#—~^y i - i : : : : : : ! V • ; • : • "• 1 i : : : : : : i 1 _; : : _ _. . . : _ : _ . . _. : _ I a. AXIAL POSITION ECS INLET FLOW AXIAL POSITION ECS INLET FLOW CO O Q. _ l < Q < DC V \ \ \ \ \ ^ r ^ 1 i I i i •• t i ' S . j . i 1 1 1 «_J1 •Jtr — B - r r - 5 — : — Z ! • / • • / • • • • / y ys ^y....:..^ r^rr... ___~p| — ^ ? - " HI jMM \ /J JUp i / '//Willi 1 '• * , 0 - - ^ A / % y / / / ' ^ > j v > 7 / ^ y y i -0^§2^ \ ^Y^ AXIAL POSITION Figure 5.12: ECS concentration field after a) 20 min, b) 40 min, c) 60 min inoculation, cin - 5 kg/m3 Chapter 5: Testing and Application of the Porous Medium Model 75 ECS INLET FLOW z o E 0) o a. _i < Q < a. AXIAL POSITION O ECS INLET FLOW AXIAL POSITION ECS INLET FLOW AXIAL POSITION Figure 5.13: ECS concentration field after a) 2 min, b) 4 min, c) 6 min inoculation, cin = 50 kg/m3 Chapter 5: Testing and Application of the Porous Medium Model 76 The ECS protein distribution at the end of the inoculation phase with each inlet concentration (i.e., c,„ = 5 or 50 kg/m3) served as a starting point for relaxation tests in which all inlet and outlet ports of the HFBR were closed. The relaxation phase can result in, after sufficient time, a reasonably uniform spatial distribution of cells and growth factors in the ECS, even if the distribution of inoculum at the end of the inoculation phase was very non- uniform. In the cases tested here, slow redistribution of protein was observed with the local concentrations reaching a relatively high degree of spatial uniformity after about 20 h or more. The process is driven by difiusion and by a weak convective flow due to the presence of osmotic pressure gradients in the ECS. The concentration field and hydrodynamics in the HFBR after 1 h and after 20 h of relaxation are displayed in Figures 5.14-5.17. The velocity vector plots, particularly those corresponding to the early stage of relaxation (Figures 5.14 and 5.16) show that the magnitude and direction of the ECS convective flow are closely associated with the magnitude and direction of the local concentration gradients. The ECS fluid travels from left to right, or towards the downstream end of the reactor, while the direction of the lumen flow remains primarily from right to left (note that the lumen velocity has no radial component and thus its distribution can be represented by a contour plot). Small areas of positive lumen velocities are located near regions of locally positive axial concentration gradients (Figures 5.14, 5.16 and 5.17). The maximum magnitude of lumen velocity correlates roughly with the concentration maximum and correspondingly decreases when the latter decreases with time. A qualitative comparison of Figures 5.15 and 5.17 indicates that inoculation with the low-concentration solution resulted in a more uniform spatial distribution of protein after 20 h of relaxation. This is also reflected by the smaller magnitude of the ECS and lumen flows in this case (Figure 5.15). It should be pointed out that, because of the no-flux condition through the boundaries, the protein transport near x = L is primarily by diffusion rather than Chapter 5: Testing and Application of the Porous Medium Model 11 1 i 1 «_ i V- | y 8 1 \ j 1 1 1 1 1 1 : 1 1 I 1 1 1 1 1: 1 I / / / / ' 1 1 1 111 1 1' I ' / / / /: : J I f / I l l M ' / 1 1 1 1 1 11111 II ' / / / / : II / / / / / / 1 • / : / J 1 1 1 III] 11 i ' 8< 6 / '• 6 - 9- / & I II uiiii i i /•••/••iii-iii- i / / / / / / / / : • /'//iii II' • / / / i i i i '• i '/////1 • i '/•/// / i / / / / / i / . / / i s / ' ? - / : 1 1 i i. AXIAL POSITION REFERENCE VECTOR: 50 * 10B m/s AXIAL POSITION AXIAL POSITION Figure 5.14: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field (MZ,-108 m/s) after 1 h relaxation following 60 min of inoculation with c,„ = 5 kg/m3. Chapter 5: Testing and Application of the Porous Medium Model 78 1 I 1 1 1 1 1 1. 1 i i i l I i i i i i J \ / ) i i i i. i i J ( ' J i î _; ^/ f • • • • • . i \ /11 y fin . y 11 V: 4 i3 2 1 AXIAL POSITION REFERENCE VECTOR: 5 0 * 1 0 " 8 m / s 0 b. ° AXIAL POSITION 1 • / / : / / / : / < 1 / 1 1 1 1 1 I \ I .^yr /. ./ 1. 7 7 i 1 ^r—t , 1—; ; / . / ; / : J • J • I X i / / I / • . / • / i I xrzxx -::X":X:"Y::K:"^:. ; yX ; 10 ; /S yS \ yf yX ;\ i , > ; ^yX ^jys • jX yX j / : >^ ; ̂ y^-^^y^ -^yX^X J^^S, | / / yr / ; >/; yX\ / / \ \\ J / / / / / \ c. AXIAL POSITION Figure 5.15: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field (uL-108 m/s) after 20 h relaxation following 60 min of inoculation with c,„ = 5 kg/m3. Chapter 5: Testing and Application of the Porous Medium Model 79 AXIAL POSITION REFERENCE VECTOR: 200*10"8m/s *^ +r ^ jr ^ > > > > > > > » » » > • ^ . AXIAL POSITION C. AXIAL POSITION Figure 5.16: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field (i/i-108 m/s) after 1 h relaxation following 6 min of inoculation with c,„ = 50 kg/m3. Chapter 5: Testing and Application of the Porous Medium Model 80 1 i I i i i "—"*~-~L^ _l """"""""••J *>>i ~— 2 1 - > i ^ 3 J \J 1 M i l l y J/I'll i* \_^^< //] ! ; ̂ y J i/\ i j ^ / j i*i / j ° ; ^^j^y/yjj ; : 4 i — 3 ^ / > ^ / : / / 1 '"̂ / / / ' ' /I ' r/ ! AXIAL POSITION REFERENCE VECTOR: 200* 10"8 m/s b. ° AXIAL POSITION Figure 5.17: a) ECS concentration field (kg/m3), b) ECS velocity field, c) lumen velocity field (wi-108 m/s) after 20 h relaxation following 6 min of inoculation with c,„ = 50 kg/m3 Chapter 5: Testing and Application of the Porous Medium Model 81 by osmotically-driven convection. Since the former transport mechanism is in this case much slower than the latter, the time for the protein distribution, for example in Figure 5.15, to reach a visually uniform state might be considerably longer than the 20 h relaxation period tested. In conclusion, inoculation with a low-concentration solution seems to be better for practical use since it can facilitate more uniform distribution of cells and growth factors over the volume of ECS. If by the time the desired average inoculum concentration is reached its distribution is not sufficiently uniform, relaxation with all ports closed will help homogenize the contents of the ECS, although this process may be fairly time-consuming and thus increase the risk of cell death due to oxygen limitations and decreases in pH. Alternatively, one might introduce the inoculum through both ECS ports or initiate the lumen flow soon after inoculation. Investigation of these options goes beyond the scope of this work but it poses no special modelling difficulty as the Porous Medium Model and the existing finite difference code are capable of handling this and even more complex cases. 5.6. Harvesting Typically, the hollow-fibre membranes in use for mammalian cell culture have low enough molecular weight cut-off values to retain the product (protein) in the extracapillary space. During harvesting, solution containing the product is collected from the ECS through its downstream port. Several different modes of harvesting are possible (Figure 5.18). In this work, two of them have been compared: the standard mode (Figure 5.18a) and the closed- lumen mode (Figure 5.18b). As mentioned before (Section 3.4), it is convenient to estimate the hydraulic conductivity, k', and protein diffusivity, D°, in the cell-packed ECS as functions of its porosity, Chapter 5: Testing and Application of the Porous Medium Model 82 Figure 5.18: Some of the possible modes of harvesting: a) standard, b) closed-lumen, c) with all ports open and equal pressures at the ECS and lumen inlets, d) with closed lumen inlet. e'ECS. Two values of e'ECS were chosen here: (i) 26%, corresponding to densely packed spheres arranged in a rhombohedral array (Bear, 1972), and (ii) 5%, an arbitrary, but more realistic, value. In each case, k* and D' were calculated from the Carman-Kozeny equation (Eq. 3.29) and from the Neale-Nader relationship (Eq. 3.30), respectively. Parameters used in the harvesting simulation study are summarized in Table 5.7. A specified concentration distribution in the ECS must be assumed as a starting point for the harvesting phase. Here, two extreme cases were considered: (i) uniform initial concentration field, c0 = 5 kg/m3, and (ii) downstream-polarized concentration field with average concentration 5 kg/m3, obtained as the steady-state solution of the one-dimensional closed-shell problem (Figure 5.19). Steady state was reached after about 19.5 h for s'ECS = 26% and 142 h for e'ECS = 5%, with the convergence criterion EPSC, defined as the maximum local concentration change with Chapter 5: Testing and Application of the Porous Medium Model 83 time (Eq. 4.4), set to 2-10"5 kg/(m3-s). One would expect similar steady-state profiles in both cases, although, at higher porosity, the concentration gradient should be steeper (since an increase in s*ECS produces larger increase in the ECS hydraulic conductivity than in the protein diffusion coefficient, see Table 5.7) and the maximum local concentration higher (because the axial ECS pressure gradients are smaller and, hence, higher osmotic pressures are required to counteract the resulting larger hydrostatic pressure differences, PL - Ps, at the downstream end). It should be pointed out that, in the 5% porosity case, the time scale for protein polarization was so large and the numerical convergence so slow that the profile shown in Figure 5.19 may not exactly correspond to steady state. However, this approximate concentration distribution is still a reasonable starting point for the harvesting simulation. Table 5.7: Summary of parameters used in the harvesting study (Gambro HFBR). Initial average protein concentration, Co Initial ECS outlet flow rate, Qs, out Initial lumen flow rate (standard mode), QL, » Membrane permeability, Lp Cell-packed ECS porosity, e'ECS ECS hydraulic conductivity, k' Protein diffusivity, D' Lumen axial permeability, &*, i Viscosity, /x 5 kg/m , (1) uniform (2) downstream-polarized 1.0000 cmVmin 600.00 cm3/min 6.4-10"15 m (1)26%, (2)5% (l)1.03-10-13m2 ( e ^ = 2 6 % ) (2)4.43-10-16m2 ( e ^ = 5%) (1) 1.9-10"11 m2/s (sECS = 26%) (2) 3.4-10-12 m2/s (e^C5 = 5%) 7.159-10"10m2 (Eq. 3.25) 0.001139Pas (water, 15°C) Chapter 5: Testing and Application of the Porous Medium Model 84 Often, in the standard ultrafiltration mode (Figure 5.18a), the ECS and lumen outlet pressures are equal (and atmospheric). However, model simulations have shown (Figure 5.20) that, because of the small hydraulic conductivity of the cell-packed ECS, this results in such low ECS flow rates that harvesting would be extremely time-consuming. Thus, it was concluded that, in order to improve the efficiency of the process, either the lumen inlet and outlet pressures should both be raised (to the same extent to maintain a constant lumen flow rate) with the ECS outlet pressure kept at one atmosphere or Ps,dn should be lowered below PL,N (e.g., by pumping the fluid out of the ECS). Similar large decreases in ECS outflow were found for the closed-lumen mode (Figure 5.18b) when changing from the higher to lower 40- 35' 30- § 25- c .2 I 20- c 8 8 15- c I «H 5H ECS IW\NIFOLD BOUNDARIES FLOW Axial position Figure 5.19: Steady-state ECS protein concentration as a function of axial position for different ECS porosities. Chapter 5: Testing and Application of the Porous Medium Model 85 porosity at a constant ECS inlet pressure. Henceforth, to allow better comparison between the different cases, all inlet and outlet pressures were set to values that would ensure the initial flow rates of 1.0000 cm3/min through the ECS and, in the standard case, 600.00 cmVmin through the lumen (Table 5.8). The ECS flow rate was usually found to decrease with time, in the most extreme case by 9% over the time period tested. This can be ascribed to the declining osmotic effects of the protein being continuously removed from the ECS. The osmotic effects are also responsible for the non-zero ECS outlet flows at zero lumen flow (Figure 5.20). 0.25 0.25 0.20 - 0.15 — 0.10 0.05 0.00 0 100 200 300 400 500 600 700 800 900 1000 Lumen inlet flow (cm/min) Figure 5.20: The ECS outlet flow as a function of the lumen flow for different ECS porosities in the standard ultrafiltration mode with Ps,dn = PL,N =1.0 atm. Chapter 5: Testing and Application of the Porous Medium Model 8 6 In all the harvesting simulation cases, the cumulative protein removal, the ECS outlet concentration, the concentration in the harvesting reservoir, as well as the ECS concentration field were determined as functions of the total fluid volume collected from the ECS (Figures 5.21 - 5.28). With a uniform initial concentration field (Figures 5.21 - 5.25), the ECS flow rate was found to change only negligibly (< 2.5%, in the most extreme case) during the harvesting period, so the relevant variables were plotted as functions of time. As can be seen in Figure 5.21, at 26% porosity almost complete protein removal is achieved after 2 h of harvesting, whereas in the 5% case, although about 10% of the protein is removed within 2 min, the maximum removal is still less than 20% after 2 h. At s'^g = 26%, there is a significant increase in the fraction of protein removed from the ECS after about 1 h in the closed-lumen case and the differences in efficiencies of both harvesting modes become more visible. Similar increases at t « 50 min in the ECS outlet and harvesting reservoir concentrations are also observed (Figure 5.22). In contrast, the standard mode curves show Table 5.8: Pressure drops corresponding to the initial ECS flow rate of 1.0000 cmVmin and (standard mode only) the initial lumen flow rate of 600.00 cm3/min. • £ECS 26% 26% 5% 5% Initial concentration field uniform polarized uniform polarized Closed-lumen Ps.in-Ps,dn (Pa) 13,779 13,168 109,697 108,767 Standard PL,O~PL.N (Pa) 4,387 4,387 4,389 4,389 PL,o-Ps,dn (Pa) 10,333 9,881 58,804 58,051 Chapter 5: Testing and Application of the Porous Medium Model 87 1. CLOSED-LUMEN HARVESTING 2. STANDARD HARVESTING f£ c s =5% 100 -90 -80 -70 -60 -50 40 H30 20 h 10 10 20 30 40 50 60 ~\—'—r- 70 80 100 110 120 Time (rrin) Figure 5.21: Fraction of protein removed from the HFBR as a function of harvesting time at different ECS porosities (uniform initial concentration field). no inflection points. To better understand the origin of the result in the closed-lumen case, it is useful to look at the transient changes in the ECS concentration field (Figures 5.23, 5.24a). As opposed to the standard mode, in which the initially uniform concentration decreases uniformly with time, a high degree of non-uniformity in protein distribution is visible here. A significant fraction (83%) of the flow entering the ECS passes into the lumina at the upstream end of the reactor and returns to the ECS downstream. Owing to protein filtering at the membrane surface, a region of maximum concentration develops upstream, shifts downstream, and eventually reaches the outlet ECS port after about 1 h of operation. This is reflected by Chapter 5: Testing and Application of the Porous Medium Model 8 8 the increased protein concentration at the outlet and the enhanced protein removal. After a subsequent 15-20 min, the outlet concentration drops to a value lower than in the corresponding standard harvesting case (Figure 5.22). However, the concentration in the outlet reservoir still remains higher than in the standard mode. In the closed-lumen case, with e'ECS = 26%, only about 17% of the fluid crosses the half- length cross-section of the reactor inside its extracapillary space. With e ^ = 5%, in both harvesting modes, the flow in the ECS is hindered to such an extent that almost all the fluid Figure 5.22: The ECS outlet concentration and the concentration in the harvesting reservoir as functions of time (uniform initial concentration field, e'ECS = 26%). Chapter 5: Testing and Application of the Porous Medium Model 89 ECS INLET FLOW ECS OUTLET FLOW z o 2 - I < Q < <r \ \ \ \ ^...i.lii. ^sS. v[o./ >y7 ^ N \ S ^ 2 / 7 / : i e s ^^^A/?i W ^ . 6 ^ ^ N^ r"i0 :' ! 1 I i i i i 1 ; 1 [ 1 4 i I i 1 I i : I " 1 1 i AXIAL POSITION AXIAL POSITION AXIAL POSITION Figure 5.23: ECS protein concentration field (kg/m3) in the closed-lumen harvesting with 'ECS 5%: a) after 5 min, b) after 30 min, c) after 60 min. Chapter 5: Testing and Application of the Porous Medium Model 90 ECS INLET FLOW ECS OUTLET FLOW AXIAL POSITION ECS INLET FLOW ECS OUTLET FLOW AXIAL POSITION ECS OUTLET FLOW * i V j . J ... i . . L.J : . ! V l :, 1 i / ! i i . ..! . . . ] ... j !. ... ;. . . j j 1 : : : ! : : i : : : : : I : : : i : : ! i : : : : : : 1 i i ! 1 ! i !"v ...u |! ll 11 ! II \ .2 . . . 1 3 iz AXIAL POSITION Figure 5.24: ECS concentration field (kg/m3) after 2 h of harvesting (uniform initial field), a) e'ECS = 26%, closed-lumen, b) e'ECS = 5%, closed-lumen, c) e'ECS = 5%, standard. Chapter 5: Testing and Application of the Porous Medium Model 91 (99.99%) travels downstream inside the fibres, except for the regions near the open ECS ports. Protein is effectively removed only from the downstream region proximal to the outlet ECS manifold, where the fluid passes back from the lumen side into the ECS. It takes just a few minutes to remove most of the protein from that region (Figure 5.25). After that time, the removal curve levels off (Figure 5.21) and the outlet concentration quickly approaches zero (Figure 5.25). In the closed-lumen case, with e'ECS = 5%, the protein in the upstream region is swept by the convective flow towards the centre of the reactor, where the concentration consequently increases. The ECS concentration contours after 2 h of harvesting are shown in Figures 5.24b and 5.24c for the closed-lumen and standard mode cases, respectively, with s'ECS = 5%. If started from a downstream-polarized concentration field, a faster and more complete protein removal can be achieved (Figures 5.26 and 5.27). After only 20 min, almost 100% and about 70% of the protein is removed at e^ = 26% and e'ECS = 5%, respectively (Figure 5.26). Plots like those shown in Figure 5.27 can be useful in estimating the time after which the protein concentration in the harvesting reservoir is above a desired level. For example, in order to obtain a 15 kg/m3 product solution, one should carry on harvesting for about 5-6 min if e^,5 = 26%, or about 20-30 s if e'ECS = 5%. The corresponding fractions of protein removed can be read from Figure 5.26 as approximately 80% and 25%, respectively. Figure 5.28 shows the ECS concentration contours as well as the ECS velocity field and the lumen velocity contours after 10 min of harvesting in the closed-lumen mode at e'ECS = 5% (polarized initial concentration field). The corresponding concentration contours after 10 min of standard harvesting look the same as in Figure 5.28a. Between 10 and 20 min, hardly any further changes in the concentration field occur. The visible high concentration ridge with a maximum above 12 kg/m3 constitutes an interesting flow division region in the closed-lumen case. Locally high osmotic pressures drive the fluid from the lumen side into this region; then, part of the fluid backflows in the ECS until it is carried by radial gradients towards the centre Chapter 5: Testing and Application of the Porous Medium Model 92 of the reactor to finally pass back into the lumina. Streamlines rather than velocity vectors are drawn in the ECS inlet and outlet regions (Figure 5.28b), where the strength of the primary hydrostatic-pressure-driven flow is two to three orders of magnitude smaller than that of the osmotically-driven flow in the central part of the ECS. The lumen velocity contours plotted in Figure 5.28c indicate that most of the fluid travels downstream near the outer wall of the cartridge (r = R). The magnitude of the maximum lumen velocity is, as expected, similar to the magnitude of the ECS inlet or outlet velocity. 5.0- 4.5- 1. ECSOlfl-LET 2. HARVESTING RESERVOIR -4 .5 -4 .0 - r - 40 —r— 50 5.0 - 1.0 - 0.5 0.0 60 Time (min) Figure 5.25: The ECS outlet concentration and the concentration in the harvesting reservoir as functions of time (both harvesting modes, uniform initial concentration field, Chapter 5: Testing and Application of the Porous Medium Model 93 100 100 Fluid volume collected from the ECS (cm3) Figure 5.26: Fraction of protein removed from the HFBR as a function of the total outflow from the ECS at different ECS porosities (both harvesting modes, polarized initial concentration field). At e'ECS = 26%, in both harvesting modes (polarized initial concentration field), the protein removal from the extracapillary space is approximately uniform (data not shown). After 2 min, the highest local concentration falls to 28 kg/m3, and after 10 min to about 9 kg/m3 (38 kg/m3 being the initially highest local value). In conclusion, the two tested harvesting modes show no significant differences in their harvesting efficiencies. Although the extent of ECS penetration by the fluid is greater in the Chapter 5: Testing and Application of the Porous Medium Model 94 closed-lumen than in the standard mode, the lumen recycle flow in the latter provides an uninterrupted diffusional supply of nutrients to and removal of metabolites from the cells. Thus, the closed-lumen mode may not be practical. At high packed cell densities, the concentration of the harvested product protein may be unacceptably low unless a sufficient degree of protein polarization in the downstream part of the ECS can be achieved between consecutive harvests. 0 2 4 6 8 10 12 14 16 18 20 Fluid volume collected from the ECS (err?) Figure 5.27: The ECS outlet concentration and the concentration in the harvesting reservoir as functions of the total outflow from the ECS at different ECS porosities (both harvesting modes, polarized initial concentration field). Chapter 5: Testing and Application of the Porous Medium Model 95 z o f= (0 2 < Q 0 ^ _ _ ^ > - ECS INLET FLOW i i : : : 1 i i i r : : i ! i ; i 01 ; . ! : I : , . . ; .: I 1 I i : : : ECS OUTLET FLOW -^\ p^ . ( I I I i ^0.1 — III ' 111 ' II \ l i 11 V l \ \ \ I \ WW \ 8 10 12 \ \ \ \ \ ____ m:.. .IL . 1\ „ V a. AXIAL POSITION REFERENCE VECTOR: 50 * 10" m/s f f I \ \ \ \ \ \ t t t / / / / / / / / • / • / • / / AXIAL POSITION ECS INLET FLOW ECS OUTLET FLOW O E W 2 - i < Q 2 l \ \ S. VJ ! ! ! „ ! ! ! is y y / 1 I \ X^4- \ \ M. —6- ; ;— 1/ / / , . . \ ! ! i ; i i ; {..../...., \ L ,..! J Z \ 1 j J . \ i ! ! ! ! i i i / A j...:::...; \ ; \ i.::.:.:..:/.: 0*1 ! I ! ! V AXIAL POSITION Figure 5.28: a) ECS concentrations (kg/m3), b) ECS velocity vectors (in the central part of the ECS) and streamlines (near the port manifolds, where the magnitude of the flow is much larger than in the central part of the ECS), c) lumen velocities (uL-l05 m/s) after 10 min of closed-lumen harvesting (s'ECS = 5%, polarized initial concn. field). Chapter 6: Conclusions and Future Work 96 Chapter 6: CONCLUSIONS AND FUTURE WORK The two-dimensional Porous Medium Model (PMM), developed here to describe the hydrodynamics and protein transport in hollow-fibre bioreactors, has been tested and applied to several situations with fundamental as well as practical implications for HFBR operation. The PMM is able to handle a variety of flow configurations, including open-shell operations, and can be relatively easily extended to include another spatial dimension or geometric design details of the HFBR. In the absence of radial terms, the Porous Medium Model reduces to the one-dimensional Krogh Cylinder Model (KCM) which was previously validated by an experimental study of protein redistribution in the closed-shell case (Taylor et al., 1994; Patkar et al., submitted; Koska, 1993a). It was found that the efficiency of the line over-relaxation procedure, employed to solve the coupled elliptic lumen and ECS pressure equations, could be significantly improved by optimizing the values of the over-relaxation parameters as well as the direction of sweep over the rows and columns of the two pressure matrices. Convergence was particularly slow for the closed-shell case and for low membrane permeabilities. In most cases, since the lumen flow was orders of magnitude greater than the ECS flow, the convergence criterion for the lumen pressure could be much less stringent than that for the ECS pressure, which reduced the computational times needed to solve the pressure equations. Since the ADI method, used to solve the time-dependent parabolic protein transport equation (coupled with the hydrodynamic equations through the osmotic pressure term), is not conservative in its standard formulation, a small time step size was required in order to minimize the resulting protein mass imbalance. Simultaneous iteration of all three equations at each time level was not necessary because the hydrodynamics did not change significantly over one time increment. Thus, the pressure solutions could be lagged behind the concentration solutions Chapter 6: Conclusions and Future Work 97 without any noticeable loss in accuracy. This was confirmed by a comparative numerical study and resulted in a several-fold reduction of computational times. Model simulations of closed-shell operations confirmed the downstream polarization of the ECS protein that occurs under dominant convective transport conditions. A two- dimensional study of this case demonstrated that, in the presence of significant radial pressure gradients in the lumen manifolds, the protein was polarized in both axial and radial directions and was distributed over a larger portion of the ECS volume than in the corresponding case where no radial gradients were imposed. When the PMM was used to study the hydrodynamics of hollow-fibre devices in the partial and full filtration modes of operation, it was found that, for membranes with permeabilities lower than approximately 10"13 m (which covers the range of most commercial ultrafiltration hollow-fibre membranes), practically all of the pressure drop between the inlet lumen and outlet ECS ports was due to the hydraulic resistance of the membrane. Although this assumption is commonly made when membrane permeabilities are determined experimentally, it breaks down at greater values of Lp. A correction plot for the Gambro hollow-fibre module geometry was developed which allows better estimates of Lp values from measured flow-rate-versus-pressure-drop data. The lumen and ECS volumetric flow rates in both filtration modes were calculated using both the PMM and the KCM (Kelsey et al., 1990). The differences in the predictions of the two models became noticeable only for Lp values above about 10'12 m and were mostly due to the fact that radial ECS flows were not included in the one-dimensional Krogh Cylinder Model. In addition, the position and area of the outflow surfaces in the KCM and PMM were different. Since the permeabilities of membranes in use for most open-shell situations of practical interest are lower than 10'12 m, the hydrodynamic predictions of the Krogh Cylinder Model should be acceptable in most cases. Chapter 6: Conclusions and Future Work 98 Simulations of the inoculation process using a Gambro HFBR with a membrane permeability of the order of 10"15 m showed that, at the end of the inoculation phase, the protein concentration distribution could be very non-uniform with most of the shell side free of protein. Using a lower-concentration inoculum solution might partially alleviate this problem. Alternatively, a relaxation phase with all ports closed could be applied following the inoculation to help homogenize the contents of the ECS by diffusion and osmotically-driven convection. However, this process might be time-consuming and a long period without lumen flow might result in oxygen starvation of the cells. It is suggested that introduction of the inoculum through both ECS ports or periodic changes of the flow direction may be more efficient ways of distributing the ECS proteins. The PMM could provide useful assistance in determining the optimum inoculation procedure. The harvesting phase of cell-packed Gambro HFBRs was also modelled. A comparison of predicted harvesting results obtained using the closed-lumen mode (both ECS ports open) and the standard mode (downstream ECS port and both lumen ports open), showed no significant differences. The rate of protein removal from the ECS and the product concentration in the harvested solution were greatly dependent on the cell-packed ECS porosity which determines the hydraulic permeability and thus the magnitude of convective transport in the shell side. Greatly increased product harvest concentrations were obtained in cases where the protein had been downstream-polarized prior to harvesting. However, recovery of high concentration harvests or complete removal of the product from the ECS of high-cell-density HFBRs might not be possible. The present version of the Porous Medium Model as well as the existing numerical code can be applied to a variety of other operating conditions and flow configurations, beside those investigated in this work. In addition, the PMM can be extended to include another dimension and concentration-induced density variations, which would allow gravitational effects to be taken into account. It can also be revise to account for the influence of the ECS manifolds on Chapter 6: Conclusions and Future Work 99 ECS protein redistribution during both normal closed-shell and harvesting open-shell operations. The ECS manifolds are essentially fibre-free and, as such, may not be treated as a porous medium described by Darcy's law. Utilization of the Navier-Stokes equations for the manifold and the Brinkman equation for the porous regions may be necessary in this case. Prediction of fluid flow and pressure distribution in the inlet and outlet lumen manifolds may also be worth pursuing, particularly if significant radial pressure gradients are created during normal closed-shell operation, or if significant fluid bypassing occurs through the lumen manifolds when one or both lumen ports are closed. Furthermore, the production phase of HFBR operation can be modelled by including a reaction rate term in the protein balance equation, where the reaction rate would depend upon the cell density. Protein leakage from the ECS into the lumina as well as transport of low-molecular-weight nutrients across the membrane may be included in respective solute mass balance equations. Variation of ECS fluid viscosity and protein difrusivity with protein concentration may also be considered. Redistribution and growth of cells is another potential study area in this field. This would require keeping track of the distributions of essential nutrients and possibly some metabolites. Also, the effect of cell density on the ECS porosity and, hence, hydraulic conductivity and protein diffusivity may be taken into account. Finally, the Porous Medium Model might be used to simulate the complete HFBR operating cycle, including inoculation, growth, steady- state production and harvesting. This could allow the optimization of operating conditions and the enhancement of reactor productivity. 100 NOMENCLATURE A total surface area of the hollow-fibre membranes (m2); dimensionless coefficient in Patankar's power-law scheme (Appendix C) Av membrane surface area per unit volume available for fluid transport (m"1) A2, A3 virial coefficients in Eqs. 3.13 and A.33 (kg'1 m3, kg"2 m6, respectively) ACCF acceleration factor in the time-marching algorithm (dimensionless) B dimensionless coefficient in Patankar's power-law scheme (Appendix C) Bj, B2, B3, B4 constants in Eqs. A.9, A. 10, B. 1 and B.2 (Pa) C, c actual concentration (kg m"3) Co, Co initial concentration (kg m"3) D diffusivity (m2 s'1) EPS convergence criterion for pressures (Pa) EPSC convergence criterion for concentration (kg m"3 s"1) / filtration fraction (Eq. A. 8) (dimensionless) fosm relationship between the osmotic pressure and concentration (Appendix A) G dimensionless parameter in the expression for krts in Table 3. lb (Eisenberg & Grodzinsky, 1988) J protein flux (kg m'2 s'1) Jo, Ji Bessel Sanctions of order 0 and 1, respectively k hydraulic conductivity, Darcy permeability (m2) K number of grid points in the radial direction Kc, Kj dimensionless hindrance factors for the convective and diffusive transport, respectively, inEq. A.36 (Koska, 1993a) L length of the ECS, permeable length of the dry fibre (m) Lp membrane permeability (m) LPiapp apparent membrane permeability (m) Mp molecular weight of protein (kgmol'1) ms salt concentration in Eq. 3.13 (mol m"3) 101 maxO n N o P Pe r ri r2 R, RHFBR R8 RL RM Rs Re Q t T TL u * u V * V V V* X *m zP function returning the maximum of the arguments in parentheses number of fibres in the HFBR cartridge number of grid points in the axial direction order of magnitude hydrostatic pressure (Pa) Peclet number (dimensionless) radial position (m) radial position of the south face of the control volume (m) radial position of the north face of the control volume (m) the HFBR cartridge radius (m) the gas law constant (J mol"1 K"1) fibre inner radius (m) fibre outer radius (m) Krogh cylinder radius (m) Reynolds number (dimensionless) volumetric flow rate (m3 s*1) time (s) absolute temperature (K) dimensionless transport modulus (Bruining, 1989) axial superficial velocity component (m/s) axial actual velocity component (m/s) radial superficial velocity component (m/s) radial actual velocity component (m/s) superficial velocity vector actual velocity vector axial position (m) axial length of the ECS manifold (m) protein charge number (dimensionless) 102 Greek letters a P AP APm Ar At Ax 8 SECS • SECS £s <i> <p Y K X A(\Pe\) V- K n p ¥ over-relaxation parameter (dimensionless) parameter in Eq. A. 17 (Taylor et al., 1994) (m"2) pressure drop (Pa) pressure drop across the membrane (Pa) increment in the radial position (m) time increment (s) increment in the axial position (m) porosity (dimensionless) fraction of the HFBR volume not occupied by the fibres (dimensionless) cell-packed ECS porosity (dimensionless) overall porosity of the HFBR (dimensionless) fluid source or sink (s'1) fraction of the HFBR volume occupied by the fibres (dimensionless) dimensionless geometric parameter (Eqs. A.7 and B.9) dimensionless membrane permeability (Eqs. A. 16 and B. 11) dimensionless parameter defined by Eqs. A. 15, A.25 and B.10 function in Patankar's power-law scheme (Appendix C) (dimensionless) fluid viscosity (Pas) 3.1415926535897... osmotic pressure (Pa) fluid density (kgm"3) sink or source of solute (kg m"3 s'1) Subscripts 0 AVG dn E U in L N out r R S up W x initial (at t = 0); lumen inlet (x = 0) average downstream east face of the control volume indices of the axial and radial positions, respectively inlet lumen lumen outlet; north face of the control volume outlet radial atr = R shell side; south face of the control volume upstream west face of the control volume axial Superscripts IT iteration counter cell-packed actual (as opposed to superficial) Other symbols (overbar) V 1 radially-averaged Nabla operator unit vector Abbreviations ADI Alternate Direction Implicit BSA bovine serum albumin ECS extracapillary space HFBR hollow-fibre bioreactor KCM Krogh Cylinder Model PMM Porous Medium Model REV representative elementary volume 105 REFERENCES Anderson D. A., J. C. Tannehill, & R. H. Pletcher, Computational fluid mechanics and heat transfer, Hemisphere Publishing Company, New York, 1984. ApelblatA., A. Katzir-Katchalsky, & A. Silberberg, "A mathematical analysis of capillary- tissue fluid exchange", Biorheology 11, 1-49 (1974). Baxter L. T., & R. K. Jain, "Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection", Microvasc. Res. 37, 77-104 (1989). Bear J., Dynamics of fluid in porous media, American Elsevier Publishing Company, New York, 1972. Breslau B. R., A. J. Testa, B. A. Milnes, & G. Medjanis, "Advances in hollow fiber ultrafiltration technology", Polym. Sci. Technol 13, 109-127(1980). Bruining W. J., "A general description of flows and pressures in hollow fiber membrane modules", Chem. Eng. Sci. 44, 1441-1447 (1989). Carman P. C , "Fluid flow through a granular bed", Trans. Inst. Chem. Eng. London 15, 150-156 (1937). Cima L. G., Anchorage-dependent mammalian cell culture in hollow fiber reactors: cell metabolism and mass transfer limitations, Ph.D. Dissertation, University of California, Berkeley, 1988. Colton C. K., B. A. Solomon, P. M. Galetti, P. D. Richardson, C. Takahashi, S. P. Naber, & W. L. Chick, "Development of novel semipermeable tubular membranes for a hybrid artificial pancreas", Polym. Sci. Technol. 13, 541-555 (1980). Drioli E., "Progress in the industrial realizations of ultrafiltration processes", Polym. Sci. Technol. 13, 291-304 (1980). Drummond J. E., & M. I. Tahir, "Laminar viscous flow through regular arrays of parallel solid cylinders", Int. J. Multiphase Flow 10, 515-540 (1984). Eisenberg S. R., & A. J. Grodzinsky, "Electrokinetic micromodel of extracellular matrix and other polyelectrolyte networks", PhysicoChemicalHydrodynamics 10, 517-539(1988). EthierC. R., "Flow through mixed fibrous porous materials", AIChEJ. 37,1227-1236 (1991). Gohl H., & P. Konstantin, "Membranes and filters for hemofiltration", Hemofiltration, Springer Verlag, 1986. Happel J., "Viscous flow relative to arrays of cylinders", AIChEJ. 5, 174-177 (1959). Harrop J. A., & J. I. T. Stenhouse, "The theoretical prediction of inertial impaction efficiencies in fibrous filters", Chem. Eng. Sci. 24,1475-1481(1969). Heath C. A., G. Belfort, B. E. Hammer, S. D. Mirer, & J. M. Pimbley, "Magnetic resonance imaging and modelling of flow in hollow-fiber bioreactors", AIChEJ 36, 547-558 (1990). Hermans J. J., "Physical aspects governing the design of hollow fiber modules", Desalination 26, 45-62 (1978). Humphrey A., S. Aiba, &N. Millis, BiomechanicalEngineering, Univ. Tokyo Press, Tokyo, 1985. Hwang T. H , & S. C. Yao, "Crossflow heat transfer in tube banks at low Reynolds number", J. Heat Transfer 108, 697-700 (1986). Jackson G W., & D. F. James, "The permeability of fibrous porous media", Can. J. Chem. Eng. 64,364-374(1986). KelseyL. J., M. R. Pillarella, & A. L. Zydney, "Theoretical analysis of convective flow profiles in a hollow-fiber membrane bioreactor", Chem. Eng. Sci. 45, 3211-3220 (1990). Kim I. H., & H. N. Chang, "Variable-volume hollow-fiber enzyme reactor with pulsatile flow", AIChEJ 29, 910-914 (1983). Kim S. S., & D. O. Cooney, "An improved theoretical model for hollow-fiber enzyme reactors", Chem. Eng. Sci. 31, 289-294 (1976). Kirsch A. A., & N. A. Fuchs, "Studies on fibrous aerosol filters - JJ. Pressure drops in systems of parallel cylinders", Ann. Occup. Hyg. 10, 23-30 (1967). Kleinstreuer C , & S. S. Agarwal, "Analysis and simulation of hollow-fiber bioreactor dynamics", Biotechnol. Bioeng. 28, 1233-1240 (1986). Knazek R. A., P. M. Gullino, P. O. Kohler, & R. L. Dedrick, "Cell culture on artificial capillaries: an approach to tissue growth in vitro", Science 178, 65-67 (1972). 107 Koska J., Protein polarization in packed hollow fibre bioreactors, M.A.Sc. Thesis, University of British Columbia, Vancouver, 1993a. Koska J., private communication, 1993b. Kozeny J., Hydraulics, Springer Verlag, Vienna, 1953. Krogh A., "The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue", J. Physiol. 52, 409-415 (1919). Kuwabara S., "The forces experienced by randomly distributed parallel circular cylinders or spheres in a viscous flow at small Reynolds numbers", J. Phys. Soc. Japan 14, 527-532 (1959). Lapidus L., & G. F. Pinder, Numerical solution of partial differential equations in science and engineering, John Wiley & Sons, Inc., New York, 1982. Libicki S. B., P. M. Salmon, & C. R. Robertson, "The effective diffusive permeability of a nonreacting solute in microbial cell aggregates", Biotechnol. Bioeng. 32, 68-85 (1988). MahonH. J., "Permeability separatory apparatus", US Patent 3,228,876, 1960. Michaels A. S., "Fifteen years of ultrafiltration: Problems and future promises of an adolescent technology", Polym. Sci. Technol. 13, 1-19 (1980). Mulder M. Basic principles of membrane technology, Kluwer Academic Publishers, Boston/London, 1991. Neale G, "Degrees of anisotropy for fluid flow and diffusion (electrical conduction) through anisotropic porous media", AIChEJ. 23, 56-62 (1977). Neale G. H., & W. K. Nader, "Prediction of transport processes within porous media: Diffusive flow processes within an homogeneous swarm of spherical particles", AIChEJ. 19, 112-119(1973). Park J. K., & H. N. Chang, "Flow distribution in the fiber lumen side of a hollow-fiber module", AIChEJ. 32, 1937-1947 (1986). Patankar S. V., Numerical heat transfer and fluid flow, Hemisphere Publishing Company, New York, 1980. Patkar A. Y., J. Koska, D. G. Taylor, B. D. Bowen, & J. M. Piret, "Protein transport in ultrafiltration hollow fiber bioreactors", in press, AIChEJ. 108 Pillarella M. R., & A. L. Zydney, "Theoretical analysis of the effect of convective flow on solute transport and insulin release in a hollow fiber bioartificial pancreas", ASMEJ. Biomech. Eng. 112, 220-228 (1990). Piret J. M., & C. L. Cooney, "Mammalian cell and protein distributions in ultrafiltration hollow fiber bioreactors", Biotechnol. Bioeng. 36, 902-910 (1990a). Piret J. M., & C. L. Cooney, "Immobilized mammalian cell cultivation in hollow fiber bioreactors", Biotechnol. Adv. 8, 763-783 (1990b). Piret J. M., & C. L. Cooney, "Model of oxygen transport limitations in hollow fiber bioreactors", Biotechnol. Bioeng. 37, 80-92 (1991). Piret J. M., D. A. Devens, & C. L. Cooney, "Nutrient and metabolite gradients in mammalian cell hollow fiber bioreactors", Can. J. Chem. Eng. 69, 421-428 (1991). Rony P. R., "Multiphase catalysis. II Hollow fiber catalysts", Biotechnol. Bioeng. 13, 431- 447 (1971). Salmon P. M., S. B. Libicki, & C. R. Robertson, "A theoretical investigation of convective transport in the hollow-fiber reactor", Chem. Eng. Comm. 66, 221-248 (1988). Sangani A. S., & A. Acrivos, "Slow flow past periodic arrays of cylinders with application to heat transfer", Int. J. Multiphase Flow 8, 193-206 (1982a). Sangani A. S., & A. Acrivos, "Slow flow through a periodic array of spheres", Int. J. Multiphase Flow 8, 343-360 (1982b). Spielman L., & S. L. Goren, "Model for predicting pressure drop and filtration efficiency in fibrous media", Environ. Sci. Technol. 2, 279-287 (1968). Starling E. H., "On the absorption of fluids from the connective tissue spaces", J. Physiol. 29,312-326(1896). Swabb E., J. Wei, & P. M. Gullion, "Diffusion and convection in normal neoplastic tissues", Cancer Res. 34,2814-2822(1974). Taylor D. G., J. M. Piret, & B. D. Bowen, "Protein polarization in isotropic membrane hollow-fiber bioreactors", AIChEJ. 40,321-333(1994). Tharakan J. P., & P. C. Chau, "Operation and pressure distribution of immobilized cell hollow fiber bioreactors", Biotechnol. Bioeng. 28, 1064-1071(1986). 109 Vilker V. L., C. K. Colton, & K. A. Smith, "The osmotic pressure of concentrated protein solutions: Effect of concentration and pH of saline solutions of bovine serum albumin", J. Coll. Interface Sci. 19, 548-565 (1981). Waterland L. R., A. S. Michaels, & C. R. Robertson, "A theoretical model for enzymatic catalysis using asymmetric hollow fiber membranes", AIChEJ. 20, 50-59 (1974). Webster I. A., & M. L. Shuler, "Mathematical models for hollow-fiber enzyme reactors", Biotechnol Bioeng. 20, 1541-1556 (1978). White F. M., Viscous fluid flow, McGraw-Hill Book Company, New York, 1991. Wolf C. F. W., "Liver tumor cells grown on hollow fiber capillaries: A prototype liver assist device", Polym. Sci. Technol. 13, 557-564 (1980). Zydney A. L., W. M. Saltzman, & K. C. Clark, "Hydraulic resistance of red cell beds in an unstirred filtration cell", Chem. Eng. Sci. 44,147-159(1989). Appendix A: Modelling Equations Based on the Krogh Cylinder Approximation 110 Appendix A: FORMULATIONS OF HFBR MODELLING EQUATIONS BASED ON THE KROGH CYLINDER APPROXIMATION Kelsey et at (1990): Hydrodynamics of open- and closed-shell, cell-free HFBRs The governing equations are: d2PL L, ^ = 1 6 ^ ( P L - P S ) (A.1) d2Ps L. 1 !*=-16*[y*>--V- (A2) subject to the boundary conditions PL=Pi.o a t x = 0 , (A.3) dP - T i = 0 a t x = 0 , (A.4) dx uL*= uLf> (1 - r2/Rl) a t x = 0 , (A.5) <N = ( l - / ) < o > (A.6) where Y = [ARAS ln(Rs/RM)+4R2sR2M-3R$ -R*M]/R*L , (A.7) PL,O is the inlet lumen pressure, u[p is the inlet centre-line actual lumen velocity, «£*0 and u*N are the radially-averaged inlet and outlet actual lumen velocities, respectively, a n d / i s referred to as the ultrafiltration fraction (fraction of the inlet flow that exits the HFBR through the ECS), i.e., f = Qs*ut/Qun- (A.8) Appendix A: Modelling Equations Based on the Krogh Cylinder Approximation 111 The solutions of Eqs. A.1 and A.2 are PL(x) =£, sinh(Xx/Z) + B2 cosh(Xx/X) + B3 x/L + B4 (A.9) Ps(x) =-5,/7sinh(Xx/Z) -Bj ycosh(kx/L) + B3 x/L + B4 (A. 10) with 4-y v>Lu?0 B2 =JB1/sinh(X)[l-cosh(X)-/(l+l/y)] (A.12) B3=Bjy (A. 13) B,=PLf)-B3 (A. 14) X = 4 V K ( 1 + V Y ) , (A.15) where K is the dimensionless membrane permeability, K = LpL2/R3L. (A. 16) Taylor et al. (1994): Radially-averaged velocities and two-dimensional ECS protein transport in a multi-fibre, closed-shell, cell-free HFBR The equations governing the quasi-steady lumen and ECS hydrodynamics are ?< , »— „ — 2Lr <y_[C(*„,x)] dC(J?M,x) — - x U l = - „ „ „ - — jig — 2 - x uL = - P « L J B - T S ^ z: (A17) U * s=i?2-V (w^~u^) ( A 1 8 ) with UL = u*,o at x = 0 and x = L, (A. 19) while the ECS protein transport equation is Appendix A: Modelling Equations Based on the Krogh Cylinder Approximation 112 dC dt = D d2C \d_( dC dx2 r d r l dr, .dC dC sdx sdr (A.20) with the initial and boundary conditions c = c0 ^ = 0 dx V ; C - Z > ^ = 0 dr a c = o dr Eq. A. 17, X = A/l6Z,/flL3(l+ V = l6Lp/Rl-l/y 1/Y) at t = 0, at x = 0 and x = L, at T = RM, at r = i?5. (A.21) (A.22) (A.23) (A.24) (A.25) (A.26) where y is defined by Eq. A. 7; C(RM,X) is the ECS protein concentration at the outer surface of the fibre; /»ffl(r,x) is the relationship between osmotic pressure and the local protein concentration (Taylor et al. used Eq. 3.13, valid for BSA); D is the protein diffiisivity. The local actual velocities are calculated as follows: 2i& u,,(r,x) = - . R 2-t\n(r/RM)-r2/Rl + l R M («L>- U L(X) ) (A.27) Vs(r,x) = (A.28) Wl ^ - r 3 + | ( i ? s 2 - ^ ) ( ^ - r 2 ) - ^ ( ^ l n ( i ? 5 / ^ ) - r 2 l n ( r / i ? M ) ) du*(x) dx In the absence of osmotic effects, Eq. A. 17 has the following solution: Appendix A: Modelling Equations Based on the Krogh Cylinder Approximation 113 u uL = LP 1+y cosh(X(X/2-x)) cosh(XL/2) + Y (A.29) Patkar et al. (in press): Radially-averaged axial velocity and protein concentration in the ECS of a closed-shell, cell-free HFBR The governing equations are 2 „ * d u — 16Z. u, 'P "LP UkLp d/^pOjdCCx) AJ X2US* RL(Rl-Rl)' n(Rl-Rl) dC dx (A.30) with the boundary conditions and u; = o d£_ d2C _ a(u's C) dt ~ dx2 ~ dx at x = 0 and x = L, (A.31) (A.32) subject to initial and boundary conditions identical to Eqs. A.21 and A.22 (with C replaced by C). A, in Eq. A.30 is defined by Eq. A.25. In this case, Patkar et al. related the osmotic pressure of protein (BSA) to its concentration through the following virial equation: f-<C) = RJlM, (C+^2C2 +^c3) (A.33) where Mp is the protein molecular weight and A2 and A3 are the virial coefficients obtained by fitting Eq. A.33 to experimental data. Koska (1993a): Coupled hydrodynamics and protein transport in a closed-shell, cell-packed HFBR In a two-dimensional formulation, the governing equations are Appendix A: Modelling Equations Based on the Krogh Cylinder Approximation 114 j 2 p jr L = 1 6 ^ [ p L( x ) - p sOUU + ns(xA,)] dx2 r 5 r R r dr d2Vs dx2 = 0 and 3t r dr |_ dr j dx [_ dxj es dx c ss dr (A.34) (A.35) (A.36) where Kj and Kc are the protein hindrance factors for diffusive and convective transport, respectively (usually assumed equal to 1); us and vs are the ECS superficial velocities in the axial and radial directions, respectively, calculated half-way between the pressure nodal points using a central difference approximation of Darcy's law; es is the packed ECS porosity; n s has the same form as/™*, in Eq. A.33. Equations A.34 and A.35 are subject to the boundary conditions PL — PL,O ¥L = PL,N dPs. dx dr 0 0 at x = 0 , at x = L, at x = 0 and x = L , at r = Rs (A.37) (A3 8) (A.39) (A.40) as well as the stipulation that the incoming and outgoing fluxes at the ECS/membrane interface be equal, i.e., 5PS LRL dr * A [pL(x)-ps(x,i?M)+ns(x,*M)] (A.41) at r = RM, where s 8 ~^rHRs/RM)-2R 2 s+R2 M (A.42) Appendix A: Modelling Equations Based on the Krogh Cylinder Approximation 115 The initial and boundary conditions for the ECS protein transport equation (A.36) are the same as those specified by Taylor et al. (1994) (Eqs. A.21 - A.24), except that Eq. A.23 is here replaced with ,3C ^ , ^ - C -K,D- dr = 0 at r = R M • (A.43) In a simplified one-dimensional (radially-averaged) version, Eqs. A.35 and A.36 become, respectively, d2P« RL Lr 7 ^ = -27^-Vrir[pL(x)-Ps(x,^) + ns(x,i?M)], a x [ii.s nM) KS with boundary conditions identical to Eq. A. 3 9, and (A.44) ac at d " dx r, ~ dC KdD—-d dx d dx T , U o — Kc — C s (A.45) subject to the same initial and boundary conditions as Eq. A.32 (Patkar et al., in press). Appendix B: Krogh Cylinder Equations with Pressure Boundary Conditions 116 Appendix B: KROGH CYLINDER EQUATIONS WITH PRESSURE BOUNDARY CONDITIONS IN THE CLOSED-SHELL, PARTIAL, AND FULL FDLTRATTON MODES The expressions for the lumen and ECS pressures in the formulation by Kelsey et al. (1990) have been re-derived here in terms of known inlet and outlet pressures, PL,O, PL,N and Psjn, rather than the filtration fraction/, to yield the following equations: P L O O =B> Ps(x) =Bt Bnh|Xjl+-4 yL + B2 cosh i)-'] + Pr L.0 sinhl X L) yL_ + B, cosh X— | — 1 7 I L + PL.O where (i) for /= 1 (full filtration mode): B, = P —P 1S,dn 1L,0 1 [X - sinh(X)]/7 + [I/7 + cosh(X)][l + cosh(X)/7]/sinh(X) l/7 + cosh(X) B2 = -B, sinh(X) (ii) for 0 < / < 1 (partial filtration mode): te.*-S..)[«>8h(X)-l] +(PL.N - i i , 0 ) [ l + cosh(X)/7] * i - B2 = [X -sinh(X)][cosh(X) - l ] / 7 + [X/7 + sinh(X)][l + cosh(X)/7] ^ - ^ 0 - 3 [ V 7 + sinh(X)] cosh(X)—1 (iii) for /= 0 (closed-shell mode): P —P rL,N rL,0 B, = X/7 + sinh(X) - [ 1 -cosh(X)]2/sinh(X) (B.l) (B.2) (B.3) (B.4) (B.5) (B.6) (B.7) Appendix B: Krogh Cylinder Equations with Pressure Boundary Conditions 117 _ 1 — cosh(X) B2 = B, sinh(X) (B.8) PL.O, PL,N and Ps,dn are the inlet lumen, outlet lumen and outlet ECS dimensional pressures, respectively. 7 is a geometrical factor, defined as follows: 7 K (A^jA^ VRM, \RMJ + 4 RM) \RM) (B.9) where RL, RU and Rs are the inner fibre radius, the outer fibre radius and the Krogh cylinder radius, respectively. X is defined as X = 4 ^ ( 1 + 1 / 7 ) where L.L2 K = -*• Rl (B.10) (B.H) is the dimensionless membrane permeability. The volumetric flow rates into and out of the hollow-fibre device are as follows: Qlout = ~ ^ Z ~ ^ (VT+cosh(X)) + B2 sinh(X)] &.«, = - : ! ^ f e ( l - c o s h ( X ) ) - J B 2 s i n h ( X ) ] (B.12) (B.13) (B.14) where n is the number of fibres, fi is the fluid viscosity and L is the reactor length. For L =0.215m, RL =1.15-10~4m, i?M =1.25-10_4m, i^=1.75-10"4m and Z, = 7.35-10~15m, we obtain: 7 = 0.677794, K = 2.2339-10~4, X =0.094062. The above parameter values correspond to the Gambro HFBR used in the body of the thesis (see Chapter 5), except that the value of Lp =7.35-10 ~15 m was obtained from the following calculation: Appendix B: Krogh Cylinder Equations with Pressure Boundary Conditions 118 1.5 m2/1.26 m2 • 6.18-10"15 m = 7.3510"15 m, where 1.26 m2 is the total surface area of the hollow-fibre membranes for the above-specified fibre dimensions, while 1.5 m2 is the surface area used in the determination of membrane permeability (see Section 5.3), which yielded Z„=6.18-10-,5m2. Appendix C: Patankar 's Power Law Scheme 119 Appendix C: PATANKAR'S POWER-LAW SCHEME In situations where either diffusive or convective transport strongly predominates, i.e. for either very small or very large absolute values of the Peclet number (Pe), the governing differential equation is usually discretized using either central differences or the upwind scheme, respectively. A more general approach proposed by Patankar (1980), which is valid for any value ofPe, is outlined below. The protein flux through each interface between two adjacent control volumes can be expressed in terms of the local concentrations in these volumes (see Figure 4.1), i.e., JE = Ui<jcE - Dx(dc/dx)E = {BEc^ - AEci+i.)Dx/Ax (C.l) Jw = «i.i.j cw ~ Dx(dc/dx)w = (BwciAi - Aw c^)DjA x (C.2) JN = rN v y cN - rNZ)r(3c/ar)N = (BN CjJ - AN cUi+l)Dr \(t-]+x + r^/A r (C.3) Js= rsv^lCs- rsDr(dc/dr)s= (BscUA- Asc,i)DrKri + r^/Ar (C.4) where Dx and Dr are the axial and radial protein diffusivities, respectively, and the subscripts E, W, N, S denote the respective faces of the (ij) cell. Each of the coefficients A and B is defined as follows: A = A (|Pe|) + max(-Pe, 0) (C.5) B = A{\Pe\) + max(Pe, 0) (C.6) where the function max() returns the largest of the arguments in the parentheses. The form of A(|Pe|) depends on the discretization scheme employed. For example, A(|Pe|)= 1 in the fully Appendix C: Patankar's Power Law Scheme 120 upwind scheme, while A(|Pe|) = 1 - £|Pe| in the central-difference scheme. Patankar has proposed the following expression, which is essential to his power-law scheme: A(\Pe\) = max(o, [l - 0.l|Pe|]5). (C.7) Equation C.7 yields values that are very close to those obtained from the exponential formula A(|Pe|)= |Pe|/[exp(|Pe|)-ll, which is exact in the one-dimensional case. However, the evaluation of A(|Pe|) from the power-law scheme is numerically more efficient since it does not include the exponential term. The Peclet numbers on the faces of the control volume are evaluated as follows: U:; AX PeE = - * L — (C.8) Uj_, : A x ^ v = -L j d— (C.9) v Ar PeN = ^ ~ (CIO) (CM) Appendix D: Source code in Fortran Appendix D: SOURCE CODE IN FORTRAN This program calculates the coupled pressure and conoaiiaalion fielde MI H I S hoMowmbre moroactor treated as a porous bed. Uniform grid, cylindrical ooordinatas are employed. Axial symmetry Is assumed, ihus Implying a 24tmenslonai oeee in which gravity forces are neglected. At each erne level, l i e Alternate Direction Implicit method ieuesd to find oonoonvatjon field, and the Hn*by4lne method is used Horaavely to find the lumen and shell pressure fields. This is s general version that allows for nofiux condition at file lumen bilevoutlet as well as for open-shell operations. IMPLICIT REAL'S (A4I.UO-Z) INTEGER HR1.HR2,MIN1.MIN2,SEC1,SEC2.HSEC1.HSEC2 INTEGER MONTH1,MONTH2,DAY1,DAY2,YEAR,YEAR.DofirV,APP,ERR INTEGER EXF,ERF,ELF,PACKED,HARV,T1MAX CHARACTERS STRIn,8TRout,STRRuO,STRAPP PARAMETER (K'IS.NMvO) PARAMETER (PACKED»0,HARV=O) PARAMETER (K0SP-1,KDISPC1,NDISP^,N0ISPC1) PARAMETER <MAXrTP=6O0O0,MAXIS=«,TIMAX=1,MAXrT«10OO0) PARAMETER (INDTR*1.INDC=1,INDPS=O.INDPOs=0,INDV=0) PARAMETER (IDATF'O.nDAT-'IO.APP*)) HARV-1 H harvesting PACKEO-I Hlhe ECS ispackedwlth colls rTDAT-tO: output Into the Z3AO.D77 file every 10 Heraaons INDTR-t IfVansientreeults to be printed INDC-1 If calculated oonoanh-ation values to be printed (those at oell centres) INDPS*1 IT pressure values to be printed (except for the oamotf o pressure) INDPOs«1 if EC8osmotlG pressure values to be printed I D A T F I If reed C,PS,PL,T,DT,ACCF,IS,CTOTAL,CTOT0,TOTFLW from a file DIMENSION PLO(0:K*1),PLN1<0:K*1) DIMENSION X(0:N+1),R(0:KVI),XDISP(0:N),RD|SP(0:K) DIMENSION C1(0:Nt1,0:K+1),C(0:Nt1,0:Kt1),0POem(1:N,1:K) DIMENSION P3(0;N*1.0:Kt1),PL(0:N*1,0:K+1) DIMENSION PS1(0:N+1,0:K>1),PL1<0:N*1,0:K+1) DIMENSION U(0:N,0:K«1),V(0:N+1>0:K>,UL(0:N,0:K>1) DIMENSION POem(0:N.1:K),CNSTR6(1:K) DIMENSION CNSTR1(1:K),CNSTTO(1:K),CNSTR2(2:K),CNSTT<4(1:K) DIMENSION ACSR(1:2,1:S),BSR<1:S,1:K),ABCLR<1:3.1:3> DIMENSION ABC8C(1:S,1:K,1:S).TS(0:MAXIS| COMMON /ALFASL/ ALFAS,ALFAS1,ALFAL,ALFAL1,EPSPS,EPSPL COMMON « P U CXPL.CNSTPL COMMON IXPSf CXPS,CXPS1,CXPS3 COMMON /CSTR1/ CNSTR1.CNSTR2 COMMON /CSTR3/ CNSTR3,CNSTR4,CXCN COMMON /VELf CNSTU.CNSTUUCNSTV.CNSTVR COMMON /PLuoV PUpO,PLupK,PLdnO,PLdhK COMMON /ORIDC/ CRCN.CR12.L.RHFBR COMMON /LOG/ LOOIND COMMON /CN3TRB CNSTRS COMMON lOSm OSM1,OSM2,OSM22A>2,A2,A3 COMMON /PECU CPew.CPne COMMON /PIBU PI COMMON /SWR/ SWR.SWR1 .SWL.SWL1 COMMON fSWCf 8WD,SWD1,SWU,SWU1 COMMON /ALT/ ALTERR,ALTERC.IRR.IRC,ISWOFF.INrn< COMMON/FI/FI.E.EL COMMON /N12/ N1.N2 COMMON /BCf BCPLO,BCPLN,BCPSup,BCPSdn COMMON /BCONC/ COup COMMON /PSupdnl PSup.PSdn COMMON /CNSTCN/ CNSTCN COMMON /CFLUXf QCIn.QCout,CTOTAL,CTOT0,CTOT1,VECS,TOTFLW COMMON /RDSP/ RDISP.R COMMON /INDR2/ INDR2 COMMON /ERROR/ ERR COMMON/EXFRU EXF.ERF.ELF COMMON /OUTFLW/ QoutUOoutS STPJn^ZaAG.DOO- STTtout^ZSAQ.DOO' 8TRR(XK23AG.R00- IF (APP.EQ.1) STRAPP^APPEND1 IF (APP.NE.1) STRAPP^-SEQUENTIAL' OPEN(UNrp-6,ACCESS=STRAPP,FILE=STRR00) CALL DOSTTM (HR1,MIN1,SEC1.HSEC1) CALL DOS0AT (MONTH1,DAY1,YEAR,DolW) DATA (TS( I ) . I I .MAXIS) > «3O.DO,12O.D0,18O.D0,2*O.D0,3O0.D0,3OO.D0r C C EXF,ERF,EIJ::ECSandlumanveloon>muHipHcalionfactorsfor C output C EXF*1000000 ERF-1000000 ELF1000 LOGIND1 .00 EPSPS*1.D-o EP8PL=1.D-6 I soouracy of PS ! aooureoy of PL EPSCDT1.D-8 taocuraoyofdC/dt ALFAS=1.asD0 ALFAS1=1.DOALFAS ALFAL1.6D0 ALFAL11.00ALFAL WRITE(B,,)'EPSP8: '.EPSPS WRITE(5,*),EPSPL: '.EPSPL WRITT^O'EPSCDT: ".EPSCDT WRITE(6,2S)ALFAS WR|TE(5,2S)ALFAL 25 FORMATf ALFA&'.FS.e) 26 FORMATf ALFAlj'.FS^) C Ph=4.D0*DATAN(1.D0) NFIBR=8123 RL=1.15D-4 RM=1.2SO-4 RHFBR=0.031500(2.DO RS-8HFBRfDSQRT(DBLE(NFIBR)) FI=RM*RM/RSIRS EPOR=1.D0J=l EL=RL*RL/RS/R8 L=0.216DO LECS°O.02400 I length of the ECS portmanifoid VECS=PI*L*(RHFBR^tHFBR4>BLE(NFIBR)*Rtir*RM) N1*LECS/L*N N2=N-N1*1 W R R T E ( S , * ) ' R A D I A L O R I D P O I N T S : ' , K W R I T E ( 5 , 7 A X I A L O R I D P O I N T S : \ N W R I T E ( 5 , * ) ' N I * . N I . ' N 2 = - , N 2 WRITE(S.*)<MAXtS=',MAXlS C 8urfAr=1.6O0 visoH.IStO-S IvisoosrtyofwateratlSdegC Lp°6.1SD-1S WRITEfVVLo^.Lp DHT-1.D-10 IF (PACKED.EQ.O) THEN DittX=0lll*EPOR DifTR=0irTEPOR/(ZD0^POR) ECSPOR*1.D0 ELSE ECSPOR=0.2600 DiftX=EPOR*2.D0^CSPOR/(3.D0-ECSPOR)*DrrT DiffR=OifrX END IF E*EPOR*ECSPOR VECS-VECS'ECSPOR WRITE(6,*)'EC8 as a fraction of ths total volumei'.EPOR WRITE(5,*)-CELL-PACKED ECS POROSITY:',ECSPOR WRITHS.'VEFFECTIVE ECS POROSITY:",E WRITE(6,*>,DIFFUSIvTTY: •,DhT WRITE<6,*rEFFECTiVE AXIAL DIFFUSIVTTY: \DHrX WRtTE(6,TEFFECTlVE RADIAL DIFFUSIVTTY: '.Dlffit C BCPLCM.DO 11 = FLUX,0 = NOFLUX«t the lumen inlet BCPLN=1.D0 11 = FLUX, 0 = NO FLUX at the lumen outlet BCPSup=0.D0 11 = FLUX. 0 = NO FLUX at the ECS upstream port BCPSdn=0.00 I I =FLUX.0 = NO FLUX at ths ECS downstream port C IF (BCPL0.Eai.D0) THEN PLupO=4672.2D0 PLupK=PLuf>0 WRITE(6,2) PLupO,PLupK 2 FORMATf PLupO^,F».2,- PLupK=\F».2) END IF IF (BCPLN.Eai.DO) THEN PLdn0=0.D0 PLOnK=PLdnO WRITE(6,3) PLdnO.PLdnK 3 FORMATf PLdnO=\F».2,- PLonK=-,F».2| END IF PSup=1087«7.2D0 C0up=1O.D0 IF(BCPSup.Eai.DO)THEN WRITEl6,4)PSup 4 FORMATf PSup^.FS.2) WRITE<5,5)C0up 6 FORMATf C0up=\F5.1) END IF PSdn=O.D0 IF (BCPSdn.EQ.1.00) THEN WRrTE(S,S) PSdn « FORMATf PSdn=\F».2) END IF C IF(IDATF.Eai)THEN OPEN(UNIT=3,FILE=STRin) I READ DATA FROM CALLINPF(IS,T,DT,ACCF,CTOTAL,CTOT0,TOTFLW,PL1,PS1,C1) CLOSE(UNrT=3| DT=0.1D0 WRITE<S,TlTERATION :\ITER WRlTE(S, ,) ,DT,s: ,,DT WRrfEIS.TACCF r'.ACCF WRtTE(S,TCavg .-.CTOTAUVECS WRITE<6.*),CTOTAL :\CTOTAL WRrTE(5,TCTOT0 :\CTOT0 WRrTE(5,VTOTFLW, m3 :-,TOTFLW ELSE PSO=0.600*(PLupO+PLdnO) Appendix D: Source code in Fortran 122 CinrWO.DO toonoantration of protein In th« HFBR at t=0 WRiTE<S,1>ClnH 1 FORMAT(lClnlt=l.F6.1> CTOTAL>ClnK*VECS I total amount of protein In HFBR CTOTO-CTOTAL I total amount of protein In HFBR at t=0 C DT11 .D0 T10W-MW.D0 DT10Ma2.D0*DT1 ACCF-(T10aU>T1)hm0M4T10M) WRiTE(5,') lOT0,«:\OT1 WRITE(S.')DT10m,. :'.DT10M WRITE(S,')-*CCF:,,ACCF END IF C MPXS-1 I MPRS-1 ! CALL PRMBS(RL.RM,RS,FI.MPXS,MPftS,PRMXS.PRMXL,PRMRS) C DX"UDBLE(N) DR*<HFBRfDBLE(K) CP«w=OX/DiffX CPna-DR/Dtfm CNSTCN«2.D0*PI*RHFBR*DX/CPn« CXCN-OrRX/DX/DX CRCN4HTR/0RIDR CPa8urfAr*LpfPI/L/RHFBRfRHFBR C R 1 2 M . D 0 - P R M R S / C P I D R CXPS>4 .DO*PRMXS7CPFOX/DX C X P 8 1 M C X P S - 1 . D 0 ) C X P S S ° - ( C X P S - 0 . 7 5 0 0 ) C X P L M . D O * P R M X L / C P I D X / D X C N 8 T P L « - 2 . D 0 * C X P L + 2 . 0 0 CNSTU-4>RMXS/vraafOX I ECS auparRddvdodty ooaffldant CNSTV>-PRMRS/viaa ! ECS auporSoial valooity ooaffldant CNSTVR~4>RMRSfviaon)R IEC8 auparfldal valodtyooaffldant CNSTULa-PRMXLMadDX I LUMEN auparfldal valodty ooaffldant C C SWR=SWR1M,8WL-SWL1^«waaptotharigTitrow-by-rowupelato C 8WL-8WL11,8WR-8WRla j : awaap to lha ton, row-by-row update C 8WR1=8WL1=1,SWR=SVVL=0:ltara»on-by-ltaratlon update C C 8W0-8WD1«1,8WU"«WU1-O: awaap downward*, eoluran-by-oolumn update C SWU=«WU1*1,SWD=3WD1-tt awaap upward^ oorumn-by-oduran update C SW014WU1>1,SWDaSWUH>: IteraBon-by-ltarabon update C C ISVrOFF-HfaoonvanjadpraaauralPLorPSIIatobaawttohad C oftlnturKiarltoratlona; C ISWOFF-0 It both praaiurai ai a Iterated unW both oonvarga; C C INITTt^tfooluninawaiplngatortoprtortorowawaapIng C INITR-llfrowawaaplngatartefirat C C INDR2«1 h* aaoond awaap ovar rowe kiatoad of awaap ovar C oohimna in oonoantration aquation C INDR2-1 INITR*0 I 8WOFF1 8WR1>O.D0 SWR*O.D0 SWL11.D0 SWL1.D0 ALTERR*) I ALTERR*1 If row awaap In alternate diraollona SWD11.D0 SWD1.D0 SWU1>O.D0 SWU-O.DO ALTERC=C I ALTERC-rlrl oolumn awaap In alternate diractiona IRR=1 I IRC ooiumn iwaapa avary IRR row twaapa IRCaO I IRR row awaapa avary IRC oolumn awaapa C I no row awaaplng If IRRaO, no oolumn awaaplng tf IRC=u) IF (IRR.EQ.O) mrm=o IF (IRC.EQ.O) IN ITR- I C Roaa=3.314O0 TEMP=28».D0 XMp>«a.DO Xma-150.D0 OSMIaRgaVTEMP/XMp OSM2*2.DCfXniaTOIp OSM22*OSM2*OSM2 Zpa-20.400 Zp2J€p°Zp A2a-&S2oTM-241D-4*Zn-3.sM0-S'Zp2 A3*2.>5O-6-1.06iD-S*Zp*1.72S0-7*Zp2 C IF(IDATF.NE.1)THEN T>O.DO DT"OT1 IS-1 TOTFLW=0.D0 END IF ITER=<P ITERD=0 C CALL GRID > (N,K,DX,0R,X,R,XDiSP.RDISP,CNSTR1,CNSTR2,CNSTR3,CNSTR4,CNSTRS) CALL PLBNDR(N,K,X,R,PL0,PLN1,PS0,PL1) IF (IDATF.NE.1) CALL INIT(N,K,X.R,CiniLC1,PS0,PS1) CALL COEFP(ACSR,BSR,ABCLR,ABCSC) IF (TS(1).EQ.0.D0) THEN CALLVELS<PS1,PL1,U,V,UL,RDISP,DX) CALL OUTPUT (PS1,PL1,U,V,UL,C1,X,R,XDISP,RDISP,0.D0,ITER,IS, > T,DT,NDISP,NDISPC,KDSP,KDISPC,INDC,INDPS,INDPOa,INDV) IS=2 END IF CLOSE (UNITaS) C 10 T*T«OT tTER=fTER+1 lTERD=fTERD+1 ERFM) WRITE(«,*),rrERATION:',fTER WRiTE(o,*)11l«E:",T CALLOSMOT(N,K,C1.POani,CPOani) CALLPRESRS(PL1>PS1.PL.PS,CPOani,ACSR,BSR,ABCLR,ABCSC.MAXITP) IF(ERR.EQ.1)THEN OPEN(UNlT=S,ACCESS*'APPEND1,FILE=STRR00) WRtTE(S,*)'MAXIMUM NUMBER OF ITERATIONS (PRESSURES)1 CLOSE(UNiT*6) END IF CALL VELS(PS,PL,U.V,UL.RDISP,DX) T0TFLWiTOTFLW4OT*Q0UtS CALL CONCN(U,V,DT,C1,C) IF (rTERD.EQ.ITD AT) THEN OPEN(UNIT-«.FILE=STRout) IWRITETHELASTOUTPUTTOAFILE I TO LATER READ DATA FROM CALL OUTF (IS,T,DT.ACCF.CTOTAL,CTOT0,TOTFLW.PL,PS.C) CLOSEtUNITM) ITERD=0 END IF CALL CLCMAX(0,N*1.0,K*1,C1,C,CMAX) WRITE*.*;*)1 CMAX/DT:-,CMAX/DT OPEN<UNrT=6,ACCESS=1APPEND,,FILE=STRR00) WR(TE(S,16HTER,T,DT,CMAX.CMAXfDT 16 FORMAT(1X,l4,,:T^,F14.2,' DT*,F14.4, > ' C M A X ' . F ^ a , 1 CMAX<DT*'.F12.6) CLOSEXQ IF (CMAX/DT.LE.EPSCDT) THEN OPEN(UNrr=6,ACCESS=-APPEND-,FILE=STRR00) WRITE<6,*)1STEADY STATE1 GOTO 200 ELSE IF (ITER.EQ.MAX1T) THEN Of>EN(UNIT*5.ACCESS=1APPEND',FILE=STRRuO) WRrTE(S,20) 20 FORMATIT MAXIMUM NUMBER OF ITERATIONS REACHED1) IF (T.GE.TS(IS).AND.IS.LE.MAXIS.AND.IS.NE.O) THEN IF(IS.LT.MAXIS)THEN IS=tS*1 ELSE IF (IS.EaMAXIS) THEN 18=4 END IF END IF 0 0 TO 200 ELSE IF (nMAX.Eai.AND.T.0E.TS(MAXIS)) THEN OPEN(UNrr=5.ACCESS=-APPEND-,FILE=STRF<O0) WRITE<5,30) 30 FORMAT^ MAXIMUM TIME REACHED1) OO TO 200 ELSE DT=0TACCF CALL SUBST(N,K,PS1,PS,C1,C) IF(INDTTt.EQ.1)THEN IF (T.GE.TS(IS).AND.IS.LE.MAXIS.AND.iS.NE.O) THEN OPEN<UNIT=6,ACCESS=-APPEND\FILE=3TRR0CI) CALL OUTPUT (PS,PL,U,V,UL,CJ(>R>XDISP,RDISP.CMAX.rrER. > tS,T,DT.NDISP,NDISPC,K0SP,KDISPC,INDC,INDPS,INDPOa,INDV) CALL FLUX (UUU,V,RDISP,DX,DT,H ARV) CLOSE<UNrr=s) OPEN(UNfTa4,FILE=STRout) IWRITETHELASTOUTPUTTO A I FILE TO LATER READ DATA FROM CALL OUTF (IS,T,DT,ACCF,CTOTAL,CTOT0,TOTFLW,PL,PS,C) CLOSEWNITM) IF (IS.LT.MAXIS) THEN 13=18+1 ELSE IF (IS.EaMAXIS) THEN 13=0 END IF END IF END IF G O T 0 1 0 END IF END IF C 200 CALL OUTPUT (PS.PL,U,V.UL,C,X,R,XDISP,RDI8P,CMAX,|TER,IS,T.DT, > N0ISP,NDISPC,KDSP,K0ISPC,INDC,INDPS,INDPOa,IN0V) CALL FLUX (UL,U,V,RDISP,DX,DT,HARV) IF(fTERD.GT.0)THEN OPEN<UNIT=4,FILE=STRout) IWRITETHELASTOUTPUTTO AFILE I TO LATER READ DATA FROM CALL OUTF (IS,T.0T.ACCF,CTOTAUCTOT0,TOTFLW,PL.PS,C) CLOSEtUNfW) END IF CALL DOSTIM (HR2,MIN2,SEC2,HSEC2) CALL DOSOAT (MONTH2,DAY2,YEARtDotW) CALL EXECTIME (HR1,HR2.MIN1,MIN2,SEC1,SEC2,HSEC1,HSEC2, > MONTH1,MONTH2,DAY1,DAY2) CLOSE(UNiT=S) STOP C END C c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . • • . . . . . . • • • C SUBROUTINE PRMBS(RL,RM,RS,FI,MPXS.MPRS,PRMXS,PRMXL,PRMRS) IMPUCrr REAL*S(A4t,L,0-Z) C CONST=RM*RMM.D0/FI FI2=FI*FI GOTO(10,20,30,40),MPXS C Happat; Tayior.Pirat.Bowan 10 PRMXS=CON3T(-DLOO(FI)-1-60(»2.D0-FI-0.600*FI2) WRITER,*)1 AXIAL PERMEABILITY: Happd; Taylor,Pl^aLBowan, GOT0100 C Drummond-Tahir.oquilatorai triangular array 20 PRMXS=CONST*(-DLOG(FI)-1^°75D0*2.D0*FI4.600*FI2) WRITE<6>*)1AXIAL PERMEABILITY: Drummond.Tahir, triangular array1 GO TO 100 Appendix D: Source code in Fortran 123 C Druramond-Tahlr.aquara array SO PRMXS>CON8T*(-0LOa(FI)-1-'n30(HZD0*FI-0.6D<rri2- > 0.06100*FIZ*TI2f<1.00*1.61»800*FI2*TI2)) WR[TE(6,*)'AXIAL PERMEABIUTY: Drunmiond.Tahlr, aquara a r ra / OOTO100 C arbrtrary valua 40 PRMXS>4.S40-10 o 40 PRMXS-f.OSD-13 WRITE(S.*)'AXI AL PERMEABIUTY: an arbitrary valua' C 100 PRMXL'0.12SD0*RL*RL*RL*RURSfRS C CON8T-«ON8T*0.600 GOTO(110,120,130,140.160,160,170),MPR8 110 PRMR8>CONST ,(-0LOO(FIHFI2-1.D0y(FI2<1.D0)) WRITE<6,*)'RADI AL PERMEABIUTY: HappaT GOTOSOO C Kuwabara 120 PRMRS=CON8T*(-0LOO(FI)-1-S00*2.D0*FI) WRITE(6,*)'RADIAL PERMEABIUTY: Kuwabara* OO TO 300 C 8anoanl-Aorlvoa,aquUataral triangular array 190 PrWRS=CONST>0LOG(FIM-4»D<H2.D0T=l-0.6O0*FI2) WRITE(S,*)'RADIAL PERMEABILITY: Sangari-Aorivoa, >Manoular array* ooToaoo C Sangani'Aorivoa.aquara array 140 PRMRS>CONST>(-DLOO(Flr.1-47»0>2.D0*FI-1.774D0*FI2« » 4.078O0T12TI) WR*TO6,YRADIAL PERMEABIUTY: SangankAorivoa, aquara a r ra / GOTOSOO C Druimond^Tahir,aqui*ataral Hangular array 160 PRMR8=CONSn-OLOO(FIH.-W?SOOt2.DCf*FI-0.5O0*F12- > 0.7M1D0*FI2*FI2) WRiTE(6,*)'RADIAL PERMEABILITY: Orummond.Tahlr,tr1anoular ar ra / GO TO 300 C Drumniond-Tahir.aquara array 100 PRMRS^:ONST*(-0LOO(Fl r-1-4763D0«<Z004.re»D0TI2y > (1.D0«0.48«2D0*FI-1.«04»D0*FI2)) WRITE<5,a)'RADI AL PERMEABIUTY: Drummond.Tahlr. aquara ar ra / GOTOSOO C arbitrary valua 170 PRMRS-3.200-10 O170 PRMRS*1.08O-13 WRtTE(5,*)-RADIAL PERMEABILITY: an arbitrary valua* C 300 WRITE<6,400) PRMXS,PRMXL,PRMRS 400 FORMATrPRMXS<',E1S.ai,'PRMXL^,E1I.W,'PRMRS^,E18.8| C RETURN END C c .................................................................... C SUBROUTINE PLBNDR(N,K,X,R.PL0IPLN1,PS0,PL) IMPUCrr REAL*8(A-H,L.O£) DIMENSION PL0(0:K>1),PLN1<0:K+1>,R<0:Kt1),X<0:N+1) DIMENSION PL(CKN*1,0:K*1) COMMON /PLucK PUipO,PLupK,PI_dnO,PLdnK COMMON IBC1 BCPLO,BCPLN,BCPSup.BCPSdn C L«X<N*1> C IF (BCPL0.EQ.1.D0) THEN PUXOWMjupO PUXK+ll-PLupK OPLup*PLupK-PLupO D O 1 0 J 1 . K PLO<J>=PLupO*fi<jyR<K*1)*DPLup 10 CONTINUE ELSE DO12J»0,K+1 PL0«Jr=O.D0 12 CONTINUE END IF IF (BCPLN.Eai.DO) THEN PLNKOf-PLdnO PLN1(K*1)-PLdnK DPLdn=PLdnK-PLdnO D O 2 0 J - I . K PLN1(JH>LdnO«R(J),R(K»1)*DPLdn 20 CONTINUE ELSE DO22J=-0,K*1 PLN1UHI.D0 22 CONTINUE END IF C IF(BCPL0.EQ.1.D0.AND.BCPLN.EQ.1.D0)THEN DOS0.M>,r">1 AoPLOfJ) PUO.Jr-A CONST-^PLNKJtnAVL D O 6 0 M . N PUI,J)=AtX(l)*CON8T 60 CONTINUE PL(N+1,JH»LN1(J) 60 CONTINUE END IF IF (BCPL0.EQ.1.D0.AND.BCPLN.EQ.0.DO)THEN DO80J=0,K+1 A*PL0<J) DO 70 M) ,N*1 PUI,Jp=A 70 CONTINUE SO CONTINUE END IF IF(BCPL0.EQ.0.DO.AND.BCPLN.EQ.1.DO)THEN D 0 1 0 0 J=«,K+1 A=PLN1(J) DO tW NO.N+1 PUW)=A • 0 CONTINUE 100 CONTINUE END IF IF (BCPL0.EaO.DO.AND.BCPLN.Ea0.DO) THEN 0O120J=0,K*1 DO110l=O,N*1 PL(l,J)-PS0 110 CONTINUE 120 CONTINUE END IF C RETURN END C C " M. . . . .M.M. . . . . . . . . . . . . . . . . . . . . . .1 . . . . . .1 SUBROUTINE GRID > (N,K.DX.DRARPCDISP.RDISP,CNSTR1,CNSTR2,CNSTR3,CNSTR4,CNSTR6) IMPUCU REAL*8(A-H,L,OZ) DIMENSION X(0:N+1),R(0:K41)JCDISP<&N).R0ISP(0:K),CNSTR6<1:K) DIMENSION CNSTR1(1:K),CNSTR2(2:K),CNSTR3(1:K),CNSTR4(1:K) COMMON /GRIDC/ CRCN.CR12.L.RHFBR COMMON /LOG/ LOGIND COMMON /VEU CNSTU.CNSTUUCNSTV.CNSTVR C X(0)=O.D0 XDISP(OH>.00 X(1H>.6D0*DX XDISPOHJX DO10l=2,N X(l)=4C(1>«OX*DBLE(l-1) XDISP(IH(DISP(1r>OX*DBLE(l-1) 10 CONTINUE X(N+1)*L C R(0H>.D0 R0I8P<0>=0.D0 R(1H>.5D0*DR RDISP(1H>R DO20J°2,K R<J>=R<1H>R*DBl_E<J-1) RDtSP(JH«>ISP(1rtOR*DBLE(J-1) 20 CONTINUE R(K*1)-RHFBR C IF (LOQIND.Eai.D0) THEN CNSTR1(1)=CR12(R(iyDLOO(R(2VR(1)) CNSTR6(1)=CH8TV(RDISP(1VDLOG(R(2yR(1)) DO30J-2.K-1 CNSTR1(J)=CR12(R(jyDLOG(R(JHyR(J» CNSTR5(J)=CNSTv7R0ISP<jyDLOG(R(J*1yR<J)) CNSTR2(J)=CR12/R<jyDLOG(R<JyR<J-1)) 30 CONTINUE CNSTR1(K)=CR12/R(KyDLOG(R(K+1VR(K)) CN8TRS<K>=CNSTv7R01SP<KyDLOG<R<K+iyR<K)> CN8TR2<K>=CR12fR<KyOLOG<R<KyR<K-1)) ELSE CNSTR1(1)aCR12'DR/R(1)*0.6D0*(R(1r»R(2)) CNSTR2(K)4R12fDRrR(K)*0.5D0*(R<K-1)4fi(K)) DO40J=2,K-1 CN3TR1(J(=CR12DR/R(J)*0.60CnR(J+1r*R(J» CNSTR2(Jp«R12IDR(R(J)'0.6D0*(R(J-1)4R(J)) 40 CONTINUE CNSTR1(K)=CR12(DRIR(K)*2.D0*R<K+1) END IF DO 60 J=1,K CNSTR3(JH}RCNa(1.DO«O.5D0*DR/R(J)) CNSTR4(J)=CRCNa(1.DO4).6O0*DR/R(J)) 60 CONTINUE C RETURN END C c .................................................................... C SUBROUTINE INIT(N,K,X,R,Cinrt,C,P80,PS) IMPUCIT REAL*8(A-H,L,0-Z) DIMENSION X(0:N+1),R<O:K*1),C(0:N+1,0:K+1) DIMENSION PS(0:N*1,0:K<f1) COMMON iN12>N1,N2 COMMON IBCI BCPLO,BCPLN,BCPSup,BCPSdn COMMON fBCONC/ COup C DO20M),N>1 DO10J=0,K*1 C(l,J)=ClnK c IF(I.LE.SO)THEN o C(l,J)=O.D0 o ELSE o C(l,J)=6.D0*Clnlt e END IF 10 CONTINUE 20 CONTINUE IF (BCPSup.Eai.DO)THEN D0 22M),N1 C(l,K+1>=C0up 22 CONTINUE END IF C DO40J=O,K*1 DO 30 l=0,N*1 P8(I.J)=PS0 30 CONTINUE 40 CONTINUE C RETURN END Appendix D: Source code in Fortran 124 c c c SUBROUTINE COEFP(ACSR,BSR.ABCLR,ABCSC) IMPLICIT REAL'8(A-H,L,0-Z) PARAMETER (K-1S.N100) DIMENSION ACSR(1:2.1:S),BSR(1:S,1:K).ABCLR(1:3,1:3) DIMENSION ABCSC(1:S,1:K,1:3) DIMENSION CNSTR1(1:K),CNSTR2<2:K) COMMON /CSTR1/ CNSTR1.CNSTR2 COMMON OCPUCXPL COMMON 0CP87 CXPS,CXPS1.CXPS3 COMMON IOCI BCPLO,BCPLN,BCPSup,BCPSdn C ACSR(1,1)-0.D0 AC8R<2.1)-CXPS-O.76O0 ACSR(1.2)*CXPS-1.D0 ACSR(2,2)*CXPS-1.D0 ACSR(1,S)-CXPSM>.7SO0 AC8R(2,S)=0.D0 C B8R(1,1)MCNSTR1(1HCXPS*3.28D0) B8R<2,1>--<CNSTR1<1>+2.00*CXPS»2.DO) BSR(3,1)*BSR(1.1) DO10J«2,K-1 BSR(1,J)-4CNSTR1<J)«CNSTR2tJ)4CXPS+S.26D0) BSR(2^)HCNSTR1(J)4CNSTR2(J,+2.D0*CXPS+2.D0) B8R<S,J)=88R<1,J) 10 CONTINUE BSR(1,K)*-<CNSTR2(KHCXPS+3.26O0) BSR(2,KXCNSTR2(K>2.D0*CXPS*2.D0) B8R(S,K)=B8R<1.K) C IF (BCPL0.Eal.D0) THEN ABCLR(1,1H>.D0 ABCLR(2,1)><9.D0*CXPL+1.D0) ABCLR(3,1)sCXPL-1.D0 ELSE IF (BCPL0.EQ.0.DO) THEN ABCLR(1,1)=0.D0 ABCLR(2,1)=-<CXPL+3.2SD0) ABCLR<3,1)=CXPL-0.76D0 END IF ABCLR<1,2)-CXPL-1.D0 ABCLR(2.2PM2.D0*CXPL+2.D0) ABCLR<3.2>>CXPL-1.D0 IF (BCPLN.Eai.DO) THEN ABCLR(1,3|"CXPL-1.D0 ABCLR(2,S)«4S.DO*CXPL»1.DO) ABCLR(3,3>=O.D0 ELSE IF (BCPLN.EQ.O.DO) THEN ABCLR<1,3>=CXPL-0.76D0 ABC LR(2,3)=HCXPLt3.2600) ABCLR(3,SH>.D0 END IF C ABCSC(1,1,1)-0.D0 ABCSC(1,1,2rMCNSTR1(1)*CXPS+3.2S0u) ABCSC(I,I,S)-CNSTRI(1) ABCSC(2,1,1)>0.D0 ABCSC(2,1,2)>-(CNSTR1(1r»2.D0*CXPS+2.D0) ABCSC(2,1,3)*CNSTR1(1) ABCSC(S,1,1)=O.D0 ABCSC(3,1,2r=-<CNSTR1<1HCXPS+3.2SD0) ABCSC(3,1,3r=CNSTR1(1) DO20J-2.K-1 ABC8C(1,J,1)=CNSTR2(J) ABC8C(1,J,2)MCN3TR1(JHCN8TR2(J)4CXPS+3.2SD0) ABCSC(1,J,3)=CNSTR1(J) ABCSC(2*I,1)*CNSTR2(J) ABC8C(2,J,2)"HCNSTR1(J>«;NSTR2<J>+2.D0-CXPS+2.D0) ABC8C(2*I,S)-CNSTR1(J) ABC8C(3,J.1)=CN3TR2«J) ABCSC(S,J,2)MCNSTR1(J)4CNSTR2(J)4CXPS*3.2800) ABCSC(3,J,S>=CN8TR1<J) 20 CONTINUE ABCSC(1,K,1)-CNSTR2(K) ABCSC(1,K,2)=-(CNSTR2<K>4CXP8+3.26D0>-BCP3up'CNSTR1(K) A8CSC(1,K,S)=0.D0 ABCSC(2,K.1)aCNSTR2(K) ABCSC(2,K,2MCN8TR2(K>+2.Du*CXPS+2.D0) ABCSC(2,K,S)*O.D0 ABCSC(S,K,1)-CNSTR2(K) ABCSC(3,K,2)=KCN8TR2<K>*CXPSt3.2S00)-BCPSdn'CNSTR1(K) ABCSC(3,K,SH>.D0 C RETURN ENO RETURN END SUBROUTINE OSMOT(N,K,C,POam,CPOani) IMPUCIT REAL*8(A4I,L,0-Z) DIMENSION C(0:N+1,0:K*1).PO*ni(0:N>1:K),CPOm(1:N>1:K) EXTERNAL FUNCTION FP C C SatupPOam C D O S O J I . K POani(0,J)-FP(C(0,J)) PO*m<N,J)=FP<C<N+1 ,J)) D O 2 0 I 1 . N - 1 PO*m<l,J)=FP(C<l,J)) 20 CONTINUE C C SatupCPOam C DO 401=1 ,N CPO«ni(l,J)=-4.D0*POam<l,J) 40 CONTINUE SO CONTINUE DOUBLE PRECISION FUNCTION FP(X) IMPLICIT REAL*8|A-H,L.O-Z) COMMON /OSMY OSM1,OSM2,OSM22A>2,A2,A3 FP=08M11D8QRT(Zp2TfT(»OSM22)-OS*l2+X+ArX,XtA3*XT<*X) RETURN END SUBROUTINE > PRESRS(PL1,P81,PL,PS>CPO«m,ACSR,BSR,ABCLR,ABCSC,MAXITP) IMPUCIT REAL*8(A-H,L,0-Z) INTEGER ERR PARAMETER (K=1S,N=100) DIMENSION PLJK0:K*1),PLN1(0:K*1),CPO*ni(1:N.1:K) DIMENSION PL(0:N+1,0:rvf1).PL1<0:N*1,0:K+1| DIMENSION PS(0:N-H,0:K+1).PS1(0:NH,0:K>1) DIMENSION AC3R(1:2,1:3),BSR(1:3,1:K),ABCLR(1:3.1:3) DIMENSION ABCSC(1:3,1:K,1:3),D(1:N,1:K) COMMON /ALFASU ALFAS,ALFAS1,ALFAt,ALFAL1,EPSPS,EPSPL COMMON (SV/RI SWR,SWR1,SWL,SWL1 COMMON/SWC/SWD,SWD1,SWU,SWU1 COMMON (ALT/ ALTERR,ALTERC,IRR,IRC,ISWOFF,INITR COMMON /N12/N1.N2 COMMON /BO BCPLO,BCPLN.BCPSup,BCPSdn COMMON /PSupdn/ PSup.PSdn COMMON/PRES IT/IT COMMON /ERROR/ ERR C C SWR=SWR1^,SWL=SWL1=tt«weep to the right, row-by-row update C 8WL=8WL1-1,SWR=SWR1=0:»wo«ptothol«ft.row-by-rowupdate C 8WR1=3WL1=1,SWR=SWL=0:rtarrton-by-ltorabon update C C 8WD=8W01=1,8WU=SWU1=O: awaap downward*, oolumn-by-oolumn update C 8WU=SWU1^,8WD=8WD1=0: »weep upward., oolumn-by-oolumri update C SWD1=SWU1=1,3WD=SWU=0: Itoration-byJtoration update C C ISWOFF=1iraoonvarg«lpr*Mura(PLorPS)iat>>b«awKohad C offmfurtharttarauona; C ia*rOFF=0rfbolhpr*eauraaaralteratadun<ilbo«ioonvarBe; C C INITR=0 If column tweaplng atari* prior to row aw—ping C INITR=1 If row aw* aping atari* llrat C INDS=0 INDL=0 IF(INITR.Eai)THEN INDSR=0 INOSCMRC ELSE INDSR=tRR INDSC-O END IF PLMAX=1.D10 PSMAX*1.D10 IT=0 C 30 IT=IT*1 WRITE(6,T PRESSURE ITERATION:',IT C C Swaap ovar rowa C IF(INDS.EaO.AND.IRR.NE.O.ANO.(INDSC.EaiRC.OR.INDSR.LT.IRR)) > THEN INDSR=INDSR*1 CALLCLCDS(N,K,PL1,CPOam,D) IF (SWR.Eai.D0.OR.SWR1.Eai.D0) THEN K1=1 K2=K KS=1 ELSE K1=K K2=1 K8=-1 ENO IF D 0 35J=K1,K2,KS CALL8WEPSR(AC3R,BSR.D,P81,PS,J,ALFA8,ALFA91) PS(0,J)=1.12SD0*PS(1,Jr4.126D0*PS(2,J) PS<Nvt,J)=1.12SD0*PS(N,JM>.125OO*PS<N-1>J) 36 CONTINUE D0 3 8 M . N 1 PS<I,0)=1.126D0*PS<I,1)-0.12S00*PS<I.2) PS(l,Kt1H1-D04CPSup)*(1.12S00*PS(l,KHI.125D0*PS(l,K-1)) > teCPSup*PSup 30 CONTINUE DOS7l=N1*1,N2-1 PS(I,0)=1.125DO*PS(I,1H>.12SOO*PS(I,2) PS<I,K+1)=1-12600-PS<I,KH>.12SD0*PS<I.K-1) 37 CONTINUE DOSSI=N2,N PS(l,O)=1.12SO(n>S(l,1H).12SD0*PS(l,2) PS(l,K41H1.D0-BCPSdn)*(1.12S00*l,S(l,KM>.126D0M>S(l,K-1)) > 4BCPSdn*PSdn 38 CONTINUE PS(0,O)=1.28ee26O0*PS(1,1)-0.140S2SD0*(PS(1>2hH>S(2,1)) > 40.01SD2SD0*PS(2,2) PS(0,Kt1H1.D0«CPSup)*(1.2o5S2SD0*PS(1,K)-0.14062SD0* > (P8(1,K-1)+PS<2,K)r«l.01662S00*PS<2,K-1)) > +BCPSup*PSup PS(N+1,OH.26662SDO*PS(N>1)-0.140e2SDO*(PS(N,2)4P8(N-1>1)) > 40.01S825D0*PS(N-1,2) P8(N+1,K+1)=<1 D0-BCP8dn)'(1.26S826D0*P8<N.K)-0.14O626O0* Appendix D: Source code in Fortran > (PS(N,K-1)«PS(N-1,K))40.01662600*PS<N-1,K-1)) » *BCPS<ki*P8dn IF (ALTERR.EQ.1) THEN IF (SWR.Eai.D0)THEN 8WR-O.D0 8WR1*).D0 8WL-1.D0 SWL11.D0 ELSE SWR»1.D0 8WR1-1.D0 8WL-0.D0 8WL1*).D0 END IF END IF IF <IND8R.EQ.IRft) IMD8C-0 CALLCLCMAX(O.N+1,0,K*1.P8,PS1,PSMAX) P8MAX»P8MAX/ALFAS WRITE(S,ar P8MAX*,PSMAX IF (PSMAX.LE.EPSPS) THEN IF (I8WOFF.EQ.1) THEN I N 0 8 - 1 ELSE IF (PLMAX.LE.EP8PL.OR.IN0L.Eai) INDS=1 END IF ELSE IF <rr.eQ.MAxm>) THEN ERR-1 IND8*1 END IF DOS0J-O.K*1 D O 4 0 M . N * 1 PS1(I^H>S(W) 40 CONTINUE SO CONTINUE ENOIF END IF C C S W M P ovar oolwnns C IF(IN08.EQ.0.AND.IRC.NE.0.AND.(IND8R.EQ.IRR.OR.INDSC.LT.IRC)) > THEN INDSC-INDSC+1 CALL CLCD8(N,K,PL1,CPO«ni,D) IF (SWD.Eai .D0.OR.SWD1.Eai .00) THEN NN1«1 NN2-N NNS-1 ELSE NN1-N NN2»1 NNS--1 END IF DO70l*NN1.NN2,NNS CALL SWEPSC(ABCSC,D,P81,P8,I,ALFAS,ALFAS1) P8(I.0)*1.12SD0*PS(I,1)-0.125D0*PS(I,2) PS(l,K*1)«1.12SOO-PS<I.KH>.126D0*PS<l,K-1) 70 CONTINUE D 0 7 1 I 1 . N 1 P8(I.K*1M1.DaOCP8up)'<1.12600T>S<l,K)-0.12S00*P8<l,K-1)) > 46CPSup*P8u|) 71 CONTINUE 0O72l>M1*1,N2-1 PS(I,K*1)»I.12SD0*PS<I,KM>.125D0*PS<I,K-1) 72 CONTINUE 0O73l»N2,N P8<l,K*1H1-D^8CP8dn)*(1.12SO(n>8(l,K)-0.12600*PS(l,K.1)) > 4SCPSdn*PSdn 7S CONTINUE D 0 7SJX),Kt1 PS(0*l)«1.125O0*PS(1,JH>.12SD0*PS(2,J) P8(N*1.JM-126O0*P8<N,JH).126O0-|>8<N-1,J) 75 CONTINUE PS(0.0)s1.286S2SOOY8(1,1M>.14002SDO*(PS(1,2HPS<2,1)) > «O.01M2SO0-pS<2,2) PS(O.K«1H1.D04CPSup)*(1.2aSS2SDO*PS(1,KH>.14082SDO* > (PS(1.K-1)<PS(2,K)HO.O1SS26D0*PS(Z,K-1)) > 4SCPSup*PSup PS(N'M,0)-1.2S5S2SDO>PS(N.1H>.140e2SOO*(PS<N.2)4P8(N-1,1)) > 4O.01SB2SD0*PS(N-1,2) P8(N*1.Kt1)-(1.DO-BCP3dn)*<1.28662SD0T>S<N,K)-0.14062SO0- > (PS(N.K-1)+PS(N-1,K))4O.O1SS2SO0*PS(N-1,K-1)) > <tBCPSiki*PSdn IF (ALTERC.EQ.1) THEN IF (SWD.Eai.DO) THEN 8 W D 4 . D 0 SWD1-0.D0 SWU=1.00 8 W U 1 1 . D 0 ELSE S W D 1 . D 0 8WD1-1.D0 8 W U 4 . D 0 8WU14>.D0 ENOIF END IF IF (INDSC.EO.IRC) INOSR=0 CALL CLCMAX(0,N»1,0,Kt1,P8,PS1,P3MAX) PSMAX»P8MAX/ALFAS WRrTE(»,T P8MAX«',P8MAX IF (PSMAX.LE.EPSPS) THEN IF( ISWOFF.Eai )THEN IN0S*1 ELSE IF (PLMAX.LE.EPSPL.OR.INDL.Eai) INDS=1 END IF ELSE IF (rr.EQ.MAxrrp) THEN ERR-1 INDS=1 END IF DO M M , K « t 0O80l=0,N»1 PS1(I,J>=PS<I,J) SO CONTINUE • 0 CONTINUE END IF END IF C IF(INDLEQ.O)THEN CALLCLCDLR(N,K,P81,PL1.CPOm,D) PU0,0)=PL1(0,0) PL(N*1,0)=PL1(N*1.0) PU0,Kf1H*L1(0,K*1) PL(N*1,K+1)»PL1<Nt1,Kt1) DO100J-1.K PL(OJHn.1(0vl) PL(Nt1,J>*PL1(Nt1,J) CALLSWEPLR(ABCLR,D.PL1,PUJ,ALFAL,ALFAL1,N.K) 100 CONTINUE D 0 1 0 1 M . N PUI,0H>.125D0*(1S.D0*PL(l,1)-10.D0>PL(l.2)+3.DO*PL(l.3)) PL(l,K+1)=O.126D0*(1S.DO*PL(l,K>-10.D0"PL(l,K-1) > 4&D0*PL(I.K-2)) 101 CONTINUE IF (BCPLO.EaO.00) THEN D0102 J«1,K PU0*l)>1.12S0O>PL(1IJ)-O.126OO*PL(2J) 102 CONTINUE PM0,OH>.12SO0^1S.D0*PL(0,1)-1O.DO*PL(0,2>fS.D0*PL(0,3)) PM0,Kt1)=0.12SO0M6.DO*PL(0,K)-10.D0TL<0.K-1) > +3.D0-PM0.K-2)) ENOIF IF (BCPLN.EaaD0)THEN DO103J=1,K PL<Nt1,J>i1.126D0*PL<N,J>-0.12SD0*PL<N-1.J> 109 CONTINUE PL<N*1,0t«O,12SO0"<16.DO*r'L<N*1,1>-10.D0*PL<N+1.2) > •3.D0TL(N*1.3)) PL(N+1,K*1)=O.126O0*(15.DO*PL(N+1,KH0.D0T>L<N*1,K-1) > •3.D0TM4NH.K-2)) END IF CAU-CLCMAX<0.N*1.O,Kt1,PL,PL1.PLMAX) PLMAX*4>LMAWALFAL WRITE(«,T PLMAX^.PLMAX IF (PLMAX.LE.EPSPL) THEN IF (ISWOFF.Eal) THEN INDLsl ELSE IF(PSMAX.LE.EPSPS.OR.INDS.Eai) INDL*1 END IF ELSE IF (IT.EQ.MAXITP) INDL=1 0O120J=0,K+1 DO110l=O.N+1 PL1(U)=PUW) 110 CONTINUE 120 CONTINUE ENOIF ENOIF C C Ch«akoonmrganoeofthepr«aaurM C IF (IND31N0L.EQ.0) OO TO 30 WRITE(e,*) IT.' PRESSURE ITERATIONS' C RETURN END C c.................................................................... C SUBROUTINE CLCMAX(N1,N2,K1,K2.A.A1,AMAX) IMPUCIT REAL*8(A-H,L,02) DIMENSION A|N1:N2.K1:K2),A1(N1:N2,K1:K2) C AMAX=0.00 0O20l=«1,N2 DO10J=K1,K2 AMAX=0MAX1(AMAX,DABS(A(I.JH>.1(I,J))) 10 CONTINUE 20 CONTINUE C RETURN ENO SUBROUTINE CLC0S(N.K.Pl,CPOm.D) IMPUCIT REAL*S(A4<,L,0-Z) DIMENSION PL(0:Nf1,0:K+1),CPO«ii(1:N,1:K),D(1:N,1:K) COMMON « C f BCPLO.BCPLN,BCPSup,BCPSdn C DO30J=1,K D<1,J)=^BCPL0"<2.D0*PL<0,J>+PL<1,J)+PL<2,J)) > -(1.D0-BCPL0)«(3.25O0T>U1.JHO.75O0-PL<2,J)) > *CPOwn(1,J> DO20N2.N-1 D(I.J)MPI-<l-1J)*2.D0T>L(I.J)*PL(l*1,J))«;PO«ii(l,J) 20 CONTINUE D(N,J)=*CPLN*(2.D0T>UN*V)+PL(N,J)+PL(N-1,J)) > -(1.D0-BCPLN)13.2SDCn>MN,J>4O.76O0*PL<N-1.J)) > •CPOwTI<N,J> 30 CONTINUE RETURN END SUBROUTINE CLC0LR(N,K,PS,PL,CPOm,D) Appendix D: Source code in Fortran IMPLICIT REAL*8(A-H,L.O-Z) DIMENSION PS(0:N*1.0:K+1),PL(0:N+1,0:K*1) OIMENSION CPO*m<1:N,1'.K).D<1:N,1:K) COMMON IXPU CXPL.CN8TPL COMMON IBa BCPL0.BCPLN.BCP8up.BCPS* C C PL, row awaap C D O S O J I . K D(1.J)-eCPL(rCN8TPL*PL(0,JH».26O0*P8(1,J)«).76O0*P8<2.J)| > -CPOani(1,J) DO20I-2.N-1 D(IJ)--<P8(l+1^H2-D0-P8<U>*P8<l-1^)>-CPOw»(l,J) 20 CONTINUE D(N,J)-BCPLN'CNSTPL*PL(N+1,J) > -<S.2600*PS<N,J>40.76D0*PS<N-1,l)>-CPO«ii<N,J) SO CONTINUE C RETURN END C c................................................. C SUBROUTINE 8WEPSR(ACSR.BSR,D0,PS1,P8,J,ALFAS,ALFAS1) IMPUCiT REAL*8<A-H,L,0-Z) PARAMETER (K-1S.N100) DIMENSION ACSR<1:2,1:S|.BSR(1:S,1:K) DIMENSION PS1(0:N+1,0:K+1).PS<0:N+1,0:K+1) DIMENSION D0(1:N,1:K),CNSTR2(2:K).CNSTR1(1:K) DIMENSION A(1:N),B(1:N),C(1:N),0SR(1:N),P<1:N) COMMON fCSTRV CN8TR1.CN8TR2 COMMON /8WPJ SWR.SWR1.SWL.SWL1 COMMON /N12N1.N2 COMMON IBa BCPLO,BCPLN,BCPSup,BCPSdn COMMON /PSupOW PSup.PSdn C C 8WR*8WR1-1,8WL-SWL1-<r. a m p I D « I « right, row-by-row update C 8WL-8WL1-1,8WR=8WR1*): awaap to t ialafL row-by-row update C 8WR1=SWL1«1.8WR»SWL=0: Itaradon-by-ttaration update C IF(J.EO,1)THEN CNST1>CNSTR1(1) A<1)»ACSR<1,1) B<1>-SSR(1,1) C<1)«ACSR(2,1) D8R(1)-O0(1.1)-SWL*CN8T1*PS(1.2)-SWR1*CNSTn>31(1,2) CONSTA>ACSR(1,2) CONSTB=8SR<2,1) CONSTC=AC8R<2,2) DO10l«2,N-1 A(l)-CONSTA BflK>ONSTB C(l)-CONSTC DSR(IH>0(I>1)-SWL*CNST1*PS<I,2)-SWR1*CNST1*PS1(I,2) 10 CONTINUE A(N)«ACSR(1,3) B<N)*8SR(3,1) C<N)-ACSR(2,3) DSR(N)-O0(N,1raWL«CNST1*PS(N.2>SWR1*CNST1*PS1(N,2) ELSE IF (J.OT.1.AN0 J.LT.K)THEN CNSTKNSTR1CJ) CN8T2=CN8TR2<J) A(1)*ACSR(1.1) B(1)-OSR(1J) C<1)*ACSR<2,1) DSR(1HJ0(1^)JWR*CNST2TS(1J-1)-SWR1*CNST1*PS1(1,J+1) > -8WL*CN8T1*PS<1,Jt1)-SWL1*CNST2*PS1(1,J-1) CONAI*ACSR(1.2) CONCI«AC8R(2.2) CONBI4SR(2*l) DO30l«2,N-1 A(l)=CONAI B < I H ; O N B I C ( I K S O N C I D 8 R ( I ) » 0 0 ( I , J > - 8 W R - C N S T 2 * P 8 < I , J - 1 ) - S W R 1 " C N S T 1 * P S 1 ( I , J * 1 ) > - 8 W L * C N 8 T 1 - P S ( I , J + 1 ) - 8 W L 1 * C N 3 T 2 * P 8 1 ( 1 , J - 1 ) 3 0 C O N T I N U E A(N)-AC8R<1,3) B(N)-BSR(^J) C(N)«ACSR<2,3) DSR(N)-O0<N,J>-SWR'CNST2*PS<N,J-1>-SWR1*CN3T1*P31<N,J*1) > -SWL^!NST1'PS(N.J*1>-SWL1*CNST2"PS1<N,J-1> ELSE IF (J.EO.K) THEN CNST1°CNSTR1(K) CNST2=CNSTR2(K) A(1)«ACSR(1,1) 8<1)-eSR(1.K)-BCPSup*CNST1 C(1)*ACSR<2,1) DSR(1)»O0(1,K)-SWR^NST2TS(1,K-1hSWL1,CNST2T>S1(1.K-1) > -BCPSup'CNST1<PSup CONSTA*ACSR{1.2) CONSTB4SR(2,K) CON8TC«ACSR<2.2) DO00l=2.N1 A(l)-CONSTA B<IH:ONSTB-BCPSUP*CNSTI C<IH;ONSTC DSR(IH>0(l,K)-SWR*CNST2*PS(l,K-1raWL1,CN3T2T>S1(I.K-1) > -BCPSup*CNST1*PSup 00 CONTINUE DOS2I-N1+1.N2-1 A(l)»CONSTA B<l)>CONSTB C(l)*CONSTC DSR(I)=00(I.K)-SWR*CNST2T>S(I,K-1)-SWL1*CNST2*PS1<I.K.1) 62 CONTINUE DO S4 M-N2.N-1 A(l)=CONSTA B(IH«ONSTB-BCPSdn*CNST1 COHSONSTC D8R(I)=00(I,K)-SWR*CNST2*PS<I,K-1)-SWL1*CNST2T>S1(I,K-1) > -BCPSdn*CNST1*PSdn M CONTINUE A(N)*ACSR(1.3) B<N>=B3R<3,K>-SCPSdn-CNST1 C(N)=ACSR(2,3) D8R<N)-O0(N,K)-SWR"CN8T2*PS(N,K-1)-8WL1*CNST2*P81(N,K-1) > -BCPSdn*CN8T1*PSdn END IF CALLTDMA(A,B,C,DSR,P,N) DO70|s1.N PS(U)=P(I)"ALFAS+P81(I.J)*ALFAS1 70 CONTINUE C RETURN END C c.................................................................... C SUBROUTINE SWEPLR(ABCLR,D0,PL1.PUJ,ALFAL,ALFAL1,N,K) IMPUCIT REAL*8(A-H.L,0-Z) PARAMETER (NN10000) DIMENSION ABCLR(1:3,1:3),D0(1:N,1:K) DIMENSION PL1(&N*1,0:K+UPL(0:N*1.0:K+1) DIMENSION A(1:NN),B(1:NN),C(1:NN),DLR(1:NN),P(1:NN) C A(1)*ABCLR(1.1} B(1)=ABCLR(2.1) C(1)-ABCLR(3,1) DLR(1)-O0(1,J) CON8TA=ABCLR(1.2) CONSTB=ABCLR(2,2) CONSTC<ABCLR<a,2) DO10l*2,N-1 A(l)=CONSTA B(l)=CON8TB C(l)=CONSTC DLR(l>=O0(l,J) 10 CONTINUE A(N>=ABCLR(1,S) B<N)=ABCLR(2,3) C(N>=ABCLR<3,S) DLR(NH>0(N,J) C CALLTDMA(A,B,C,DLR,P,N) D O 7 0 M . N PMW)»P(I)*ALFAL4PL1(W)'ALFAL1 70 CONTINUE C RETURN EN0 C c..............>..................................................... C SUBROUTINE SWEPSC(ABCSC,D,PS1,PS.I,ALFAS.ALFAS1) IMPUCIT REAL*8(A-H,L,aZ) PARAMETER (K=18,N=100) DIMENSION ABCSC(1:3,1:K,1:3),D(1:N.1:K) DIMENSION P81(0:N+1,0;K+1).P8(0:N*1.0:K*1) DIMENSION A(1:K),B(1:K),C(1:K),D8C(1:K),P(1:K) DIMENSION CNSTR1(1:K),CNSTR2(2:K) COMMON/CSTRIf CNSTR1 ,CNSTR2 COMMON IXW CXPS.CXPS1.CXPS3 COMMON /8WC/ SWD,SWD1,SWU,SWU1 COMMON IBa BCPL0,BCPLN,BCP8up,BCP8dn COMMON /PSupoW PSup.PSdn COMMON /N127N1,N2 C IF ( I .Ea i )THEN D O 1 0 J 1 . K A(J)=ABCSC(1,J.1) B(J>=ABCSC(1,J.2) C(J)=ABC8C(1,J.3) DSC(JK)(1^)4SWU*CXPS3*PS(2,J)+SWD1*CXPS3*PS1(2,J) 10 CONTINUE DSC(K)=0SC(K>-BCP8up*CNSTR1(K)*PSup ELSE IF (I.0T.1.AN0.I.LE.N1) THEN DOS0J=1,K A<J>=ABCSC<2A1> B«J)=ABCSC(2^,2) C(J)=ABC8C(2,J,3) DSC<J)<Mt,Jr>«WD,CXPS1*PS<M,Jr*8WD1"CXPS1*PS1(l*1,J) > «WU-CXPS1*PS<H1J>46WU1*CXPS1«PS1(I-1^) 30 CONTINUE B(K)=8<K>-BCPSup*CNSTR1<K) DSC(K)=OSC(K)-BCPSup*CNSTR1(K)*PSup ELSE IF (I.GT.N1.AND.I.LT.N2) THEN DO40J*1,K A<J>=ABCSC(2,J.1) B(J>=ABCSC(W,2) C(JH=ABC8C(2.J,3) 0SC<JH>(I,J>4SWD*CXPS1*PS<I-1,J)+SWD1'CXPS1*PS1(I+1,J) » <6WU'CXPS1*PS(lt1,Jr»SWU1*CXPS1«PS1<l-1,J) 40 CONTINUE ELSE IF (I.QE.N2.AND.I.LT.N) THEN DO60J=1,K A(J)=ABCSC(ZJ,1) B(J)=ABCSC(2.J,2) C(J)=ABCSC(2,J,3) DSC(J)=0(l,J)+3W0"CXP31*PS(l-1,J)+SWD1,CXP31'PS1(lt1,J) > 4SWU*CXPS1*PS(l<f1,J>«eWU1*CXPS1*PS1(l-1,J) 60 CONTINUE B(K)=S(K)-BCPSdll*CNSTR1(K) DSC(KH>SC(K)-BCPSdn*CNSTR1(K)«PSdn ELSE IF (I.EQ.N) THEN DOS0J=1,K A(J)sABCSC(3*l,1) B(J)=ABCSC(3,J,2) C(J)=ABCSC(3,J,3) DSC(J)=0(N,J)*8WD'CXP83*PS(N-1.J>+SWU1*CXPS3*PS1(N-1,J) Appendix D: Source code in Fortran SO CONTINUE DSC(Kp«SC(K)-BCP8<ki*CNSTR1(K)>PSdn END IF CAU_TDMA(A,B.C,DSC,P.K) D O 7 0 J 1 . K PS(I,J>»P<J)'ALFAS+P81(I,J)'ALFAS1 70 CONTINUE C RETURN END C c..^.....— C SUBROUTINE VEL8(PS,PL.U.V,UL,RDI8P,DX) C C CalauMMECSandkOTMnavarllcialiMlodliaa C 0nttt«oait irsofl i«o«M«) C IMPLICIT REAL*S(A-H.L,0-Z) PARAMETER <K-18,N=100> DIMENSION P8(0:N*1,0:K+1),U(0;N,0:K+1),V(0:N+1.0:K) DIMENSION PU0:NH.O:K+1).UM0;N,0:Kt1),CN3TR6<1:K) DIMENSION RDISP(0:K) COMMON IW&J CNSTU,CNSTUL,CN8TV,CNSTVR COMMON /LOO/ LOOIND COMMON /CNSTRST CNSTR5 COMMON IBCI BCPL0,BCPLN.BCP8up,BCPSdn COMMON /PSupcW P8ufi.PSdn COMMON/N12/N1.N2 COMMON/PIBL/PI COMMON /OUTFLW/ QoutL,QoutS C V(0,0H>D0 V(Nt1,0H>.D0 DOSJ*1,K-1 Vt0^IHLOOIN0*CN8TRS(JH<1.004.OOIND)*CNSTVR)* > (PS(0,J+1>-PS(0,J» V<N*1,J)«<LOOIND>CNSTR6(JW1.DCM.OGIND)*CN8TVR)* > (P8<N+1.Jt1)-PS<N+1,J» S CONTINUE V(0,K)»eCP8up*(LOOIND*CNSTR6(KH<1.D0-LOGIND)*2.00-CNSTVR)* > (PSup-PS(0,K» V(N<f1,K)=SCPSdri*(LOGIND*CNSTRS(KH<1.D04jOGIND)<2.D0*CNSTVR)< > (PSdn-PS(N+1,K)) C DO 20 M . N V<l,0H>.D0 DO10J-1.K-1 V(UHLOOIND*CN8TR6(J)4<1.DIW.OOIND)"CN8TVR)* > (PS(W*1(-P8<W)) 10 CONTINUE V(l,KH).D0 20 CONTINUE IF (BCPSup.Eai.D0) THEN D 0 2 1 I 1 . N 1 V(l,KHLOGIND*CNSTR5(KW1-D04.OGIND)*2.D0*CNSTVR)* > (PS(I,K+1)-PS(I.K)) 21 CONTINUE END IF IF (BCPSdn.EQ.1.D0) THEN D O » l * N 2 , N V(l,K)a<LOGIND*CN8TR6<KH<1.004.OQIN0)*2.D0*CNSTVR)* > (PS(I,K*1)-PS(I,K)) 21 CONTINUE END IF C C D 0 4 0 J 1 . K U(0,JHJ.D0 UL<0.J)=eCPLO*2.DO"CNSTUL*(PL<1,J>-PL(0,J)) 0 O S 0 M . N - 1 U(l,J)-CN8TU'(P8<lt1,J)-P8<l,J» UL(I.J)=CNSTUL,(PL(I*1.J)-PI-<I,J)) SO CONTINUE U<N,JM).D0 UL(N,J)=8CPLN*2.D0*CNSTUL*(PL(N+1,J)-PL<N,J)) 40 CONTINUE C U<O.0H>.DO U(0.K+1)>O.D0 UL(0,0)>BCPLO*2.DO*CNSTUL*(PL(1,0)-PL(0,0)) UU0,K*1)-SCPUr2.D0*CNS'rUL*(PL(1,K+1)-PL(0,K*1)) D 0 2SI=1,N-1 U(I,0)>CN8TU*(PS(I+1,0)-PS(l.0)) U(I.K+1>-CNSTU'(PS<l+1,KH)-PS<l,Kt1)) UL(I.O)-CNSTUL*(PL(K1,0H>L(l,O)) UU.I,K+1)-CNSTUL*(PL(I*1,K*1)-PL(I,K*1)) 25 CONTINUE U(N,0)>0.D0 U(N,K+1H).D0 IF (BCPSup.Eai.D0) THEN IF(0ABS(U(N1.K+1)).aT.1.D2'DABS(U<N1+1,K»1)))THEN U(N1,K+1)a0.5D0*(U(N1-1,K+1)«j(N1<f1,K+1)) END IF END IF IF (BCP8dn.Eai.D0) THEN IF(DABS(U(N2-1,K*1)).GT.1.D2*0ABS(U(N2-2,K+1)))THEN U(N2-1,K*1)-O.500*<U<N2-2,K*1)4U<N2,Kt1)) END IF END IF UL(N,0)*BCPLN*2.00*CN8TUL*(PL(N+1,0H>L<N,0)) UL(N,K»1)-8CPLN*2.D0*CNSTUL*(PL<N<l<1,K*1)-PL<N,K<>1)) C QoutM.DO IF(BCPSdn.EQ.1.D0)THEN CONST2=RDISP<K)*2.DO*PI*DX DO80l*N2-1,N Qout8=Oout34CON8T2*V<I.K) 80 CONTINUE END IF RETURN END SUBROUTINE TDMA(A,B.C,D,X,N) IMPUCIT REAL'S) A-H, l_Oi ) PARAMETER (NN10000) DIMENSION A(1:N),B(1:N),C(1:N),D(1:N),X(1:N),P(1:NN),Q(1:NN) C p (D^ ; ( i yB( i ) Q<iH>(iyB<i) DO10l*2,N DEN=A(I)*P(I-1)+B(I) P0F>-C(iyDEN Q(l)=<D<l)-A(l)*Q(l-1))fDEN 10 CONTINUE X(N)=Q(N) DO20NN-1.1.-1 X(IH-(I)*X(I+1HQ(I) 20 CONTINUE C RETURN EN0 C SUBROUTINE CONCN(U,V,DT.C1,C) IMPUCIT REAL*8(A-H,L,0-Z) PARAMETER (K=18,N=100) DIMENSION C1(0:N+1,0:K+1),C(0:N+1,0:K+1),C2(0:Nt1,0:K+1) DIMENSION AAN(1:N),BBN(1:N),CCN(1:N),DDN(1:N),CN(1:N) DIMENSION AAK(1:K),BBK(1:K),CCK(1:K),DDK(1:K),CK(1:K) DIMENSION U(0:N,0:K+1),V(0:N+1.0:K>,CNSTR3<1:K),CNSTR4<1:K) DIMENSION A«<1:N,1.K),B«<1:N,1:K),Aw<1:N.1:K),Bw<1:N,1:K) DIMENSION An(1:N,1:K).Bn(1:N,1:K),A»(1:N.1:K),Ba(1:N,1:K) DIMENSION RDISP(0:K),R(0:K+1),RDISP(0:K) COMMON /CSTRS/ CNSTR3,CNSTR4,CXCN COMMON /PECU CPaw.CPlM COMMON/FI/FI.E.EL COMMON /N12fN1,N2 COMMON IBa BCPL0,BCPLN,BCP8up.BCPSdn COMMON /BCONC/ COup COMMON /CNSTCN/ CNSTCN COMMON /CFLUXf QCIn.QCoutCTOTAL,CTOT0,CTOT1,VECS,TOTFLW COMMON /ROSPf RDISP.R COMMON /INDB2/ INDR2 C C2DT"€^.DBfDT DO20J=1,K JMFJ-1 DO10l*1,N IMM-1 IF(I.NE.N)THEN P*cE?CP<w*U(l,l) CNST=1.DO-0.1DO*DABS(PaoE) AP«=OMAX1<0.D0,CN8T*CN3T*CN8T-CNST>CN8T) Aa(l,J)=>AP««OMAX1(-P«cE,0.D0) B4M)*AP«*DMAX1(P«cE,0.D0) ELSE A«N,JH>.D0 Ba(N,J)-0.D0 END IF IF(I.NE.1)THEN PortV=CP«WU<IM,J) CNST=1.OO4.1D0*OABS(PacW) APw*OMAX1(0.D0,CNST*CNST*CNST>CNST*CNST) Aw<l,J)=APw«>MAX1<-P«=W,0.D0) Bw(l,J)=APwH>MAX1<PeoW,0.D0) ELSE Aw<1,J>=O.D0 Bw(1.J)=0.D0 END IF IF(J.NE.K.OR.(J.EQ.K.AND.BCPSup.EQ.1.D0.AND.I.LE.N1). > OR.(J.EQ.K.AND.BCPSdn.EQ.1.D0.AND.I.OE.N2)) THEN P.oN=CPn.*V(l,J) IF (J.EQ.K) P«oN=0.5O0«PecN CNST=1.D04.1DO*DABS(P«cN) APn4MAX1(0.D0,CNST*CNST*CNST*CNST*CNST) An(l,J)=APn+OMAX1<-P«cN,0.D0) Bn(l,J)**Piv»OMAX1(P«cN,0.D0) An(l,K)=2,D0*An(l,K) Bn(I.K)°2.DO*Bn(l,K) ELSE An(l,K)=0.D0 Bn(l,K)=O.D0 END IF IF(J.NE.1)THEN P»cS=CPn.*V(l,JM) CN8T=1.004.1D0*DABS<PeoS) APa=OMAX1(O.D0,CNST*CNST*CN8T*CN8T*CNST) Aa<l,J)«AP»«OMAX1(-PaoS,0.D0) B^Up=APWDMAX1(P«sS,0.DO) ELSE Aa(l,1)4>.D0 B4l,1)=O.D0 END IF 10 CONTINUE 20 CONTINUE C D 0 1 1 l=0,N*1 D013J=0,K*1 C2XI,J)=C1(U) 13 CONTINUE 11 CONTINUE C C r C DO100J=K,1,-1 Appendix D: Source code in Fortran 128 CON4»CNSTR4(J) CON3-CNSTRS(J) D O 6 0 I 1 . N AAN(l)»-CXCN"Bw<l,J) BBN(l|-C20T4CXCN*(Aw(l^>«S«(U)KA^U)*CON*M!n(U)*CON3 CCN(I)--CXCN*A«<I,J) D0N(l)>ea(U)*CON4*C2(I.J-1HC 20T"C2<I,J) > •An(l,J),CON3*C2<l,Jt1> 60 CONTINUE CALLTDMA(AAN,BBN,CCN,DON,CN,N) DO 70 M . N C2(U)-CN(I) 70 CONTINUE CJ(0^IH).12SD0^».DO*C2(1,J)-C2(2,J)) C2<N<-1,J)-0.12600"<».D0-C2<N,J)-C2<N-1,J)) 100 CONTINUE DO110P-NTM.N C2(l,OH>.12SDO^a.DO*C2(l,1M:2(l,2)) C2(l,K*1)-0.126OO*(*.DO*C2(l,K>C2{l,K-1)) 110 CONTINUE D O H 1 M . N 1 C2(l,0)-0.125D0*<e.D0*C2<l,1>-C2<l,2» C2(l,K«1)-(1.DO«CPSiip)*al2SD0^9.DO*C2(l,KH:2(l,K-1)) > «ecpsup*coup 111 CONTINUE C2(0,0)4>.016a26DO*(S1.DO*C2(1,1)4.DO*(C2(1,2>4C 2(2,1))* > C2da» C2(0,K«1H1-D04CPSup)*0.016S2SD0*(S1.DVC2(1,K) > •*.D0*(C2(1,K-1)+C2(Z.K))4C2(2>K-1)) > 4SCPSup*C0up C2(N+1,K*1H>.015B2SOO*<81.DO*C2<N,K> > -«.D0"<C2(N,K-1>4C2<N-1,K))*C2<N-1,K-1)) C2(N+1.0)^016S26DO*(S1.DO*C2(N,1)4.DO*(C2(N,2)4C2(N-1,1))<» > C2(N-1,2» C 0O132MO.N+1 D0133.M>,K+1 C<I,J)-C2<I,J) 1SS CONTINUE 132 CONTINUE C C 2row«WMp C IF ( IN0R2.Ea i )THEN C DO400J-K.1.-1 CON4-CNSTR4<J) CON3*CNSTRJ(J) DO 4601-1,N AAN(l )^0(CN*Bw(U) BBN<l)^2DT*CX0N-<Av«<l,JHB^I,J))+A^I.J)-CON4*fin(I.J)"CON3 CCN(I>>-CXCN*A<(I,J) DDN(l)*«<I.J)*CON4"C<l,J-1>4C20T*C<l,J) > •An<llJ)"CON3"C<I.J*1) 460 CONTINUE CALLTOMA(AAN,BBN,CCN,DDN,CN,N) D O 4 7 0 M . N C ( U H ! N ( I ) 470 CONTINUE C(0,J)=O.126D0*<».D0"C<1,J)-C<2,J)) C(N+1,JH).126O0-(».DOM(N,J)-C(N-1,J)| 400 CONTINUE DO410l«N1*1,N C(I,0H).126D0*(>.D0*C(I,1H:(I.2)) C(I.K+1)-O.12SOO*(».D0-C(l,K)-C(l,K-1)) 410 CONTINUE 0 O 4 1 1 M . N 1 C(I,0)^I.126D0^(.D0*C(I,1)-C(I,2)) C(l,K*1)-<1.DO-BCPSup)"0.12SOO*<».D0-C<l,K>-C(I,K-1» > +BCPSup*C0up 411 CONTINUE C ELSE C C oohimn«wMp C D O 2 0 0 I 1 . N D O 1 6 0 J 1 , K CON4=CN8TR4<J) CON3>CNSTR3(J) AAK(J)»-e^U)*CON4 BBK(J)-C2DT*CXCN"(A~(l,J>+e^l,J)(+A«<I.J)*CON4*Bn(I.J)'CON3 CCK(J)«-An(l,J)*CON3 DDK(JH«XCN*(Bw<l,J)-C(l-1,J)*A^I,J)*C(l*1,J» > +C2DT*C(I,J) 160 CONTINUE CCK(KH>.D0 DDK(K>ODK(K)+An(l,K)*CON3*C(l,K*1) CALLTDMA(AAK,BBK,CCK,DDK,CK,K) D O 1 7 0 J 1 . K C(WH3K«J) 170 CONTINUE C(I,0)M>.126DO*(«.DO*C(I,1)-C(I>2M C(I.K*1H>.126DO*(S.D0*C(l,K)-C(I.K-1» IF (LLE.N1.AND.BCP8up.Eai.D0) C(l,K+1H;0up 200 CONTINUE DO 210 J=1.K C<0,J)»O.12SD0"<».D0*C<1,J>-C<2.J» C(N*1^)>4.126D0*(9.D0*C(N,JH:(N-1J» 210 CONTINUE C END IF C C(0,OH>.01S826DO*(81.00*C(1,1)4.DO*(C(1,2)4C(2,1)>» > C(2,2» C(0,K*1H1.004CP8up)*0.015S2500*(81.00*C<1.K) > 4.00*(C(1.K-1)+C(2.K))+C(2,K-1)) > «BCPSup*COup C(N+1,0H>.016826D0*(81.D0>C(N,1)4.D0*(C(N,2)4C(N-1,1)>* > C(N-1.2)) C(NH,K+1H>.016S26DO*<81.D0*C<N,K) > 4.D0*(C(N,K-1>«C(N-1,K)HC(N-1,K-1)) C IF (BCPSup.Eai.D0)THEN QCIn-O.D0 DO300l*1,N1 QCIn=QCin-Bn<l,K)*C(l,K)+An(l,K)aC(l,K+1) $00 CONTINUE QCIn=CNSTCN*QCin CTOTAL=CTOTAL*DT*QCIn END IF C IF (BCPSdn.Eai.OO) THEN QC«n>O.D0 0O310l=H2,N QCoul=OCouHBn(l,K)*C(l,K)Jtfi(I.K)*C(l,K+1) « I 0 CONTINUE QCotrt=CN8TCN*QCout CTOTAL=CTOTAL-OT"aCo<lt END IF C CTOT1-0.D0 D O 3 0 4 . M . K CONST=(RDISP(J)*RDISP{JMH>ISP(J-1)>RDISP(J-1)) > «RDISP(KyRDISP(KyDBLE(N) DO 306 M . N CTOT1=CTOT1+C(l,J)*CONST 306 CONTINUE 304 CONTINUE WRITE(«,*C CTOT1 =f.CTOT1 C IF (BCPSdn.Eai.D0.OR.BCPSup.Eai.D0) THEN WRITER*)* AVERAGE CONCN =*,CTOTAL/VECS END IF C RETURN END C c....................................... C SUBROUTINE FLUX (UL.U,V,RDISP,DX,DT,HARV) IMPLICIT REAL*8(A^.L,0-Z) INTEGER HARV PARAMETER (K=18,N=100) DIMENSION U(0:N,0:K*1),V<0:N+1,0:K),UL(U:N,0:K*1) DIMENSION RDISP(0:K) COMMON fN12IN1.N2 COMMON (PIBU PI COMMON IBCI BCPLO,BCPLN,BCPSup,BCPSdn COMMON /CFLUX/QCta,QCoutCTOTAL.CTOT0,CTOT1,VECS,TOTFLW COMMON (OUTFLVW QootUQootS C WRITER*)" WRITE(«,")" C QlnL-O.DO OinS=0.D0 C IF (BCPL0.Eai.D0) THEN DO10J=1,K OjnLOiaL4PI*UL(0,J)* > (R0ISP(J)*RDISP(JH<DISP(J-1)*RDISP(J-1)) 10 CONTINUE WR(TE(6,101) OJnUQinL"6.D7 WRITE«,101) OJnL,OJnL*S.D7 101 FORMATf INLET LUMEN FLOW. m3/«^,D24.18, > ' r.F10.6,'mUr«in)-) END IF C IF (BCPSup.Eai.D0) THEN CONST2=RDISP(K)*2.D0*PI*DX D O 2 0 M . N 1 OJnS^3inS4:ONST2*V(l,K) 20 CONTINUE WRITE<6,102) OJnS,QlnS*e.D7 WR[TE(S,102) OinS,OjnS"8.D7 102 FORMATf INLET ECS FLOW, m3/i--,D24.18. > ' f,F10.6,'mL/mlnn END IF C IF (BCPSdn.Eai.D0) THEN WRITE(6,103) QoulS,QoutS-8.D7 WRITE(e,103) Qout3,QoutS-0.D7 103 FORMATf OUTLET ECS FLOW. ra3/«=',D24.18, > ' r,F10.6.-|«L/niinn END IF C IF (BCPLN.Eai.D0) THEN Qoutt_=«.D0 DOMJ=1.K QotltL=OoutL*PI*UL(N.J)* > (RDISP(J)>RDISP(J)-RDISP(J-1)*RDISP(J-1)) • 0 CONTINUE WRITE(S,104) Qou1L,QoutL*8.D7 WRITE(<,1M) QoutL,QoutL*8.D7 104 FORMATf OUTLET LUMEN FLOW, m3l«=',D24.1S, > ' r,F10.6.-mL/min» END IF C WRfTE(5,a)'TOTAL FLUID PASSED THROUGH ECS, mL:\TOTFLW1.D8 WRITE(6,*)TOTAL FLUID PASSED THROUGH ECS, mU'.TOTFLWI.DS C Nh«lf=*/2 QLh«lf=0.D0 DO60J=1,K QUl«lf=QLh«H*P|-UL(Nh«H.J)* > (RDISP(J)*RDISP(J)-RDISP(J-1CRDISP(J-1)) 50 CONTINUE IF(BCPSup.Eai.D0)THEN WRITE(6.106) QLhairS.D7.QLhrfffOjnS1.D2 Appendix D: Source code in Fortran WRtTE(a,106) QLhdre.D7,QLhaHtOJnS*1.D2 106 FORMATf HALF-LENGTH LUMEN FLOW, mL/™ln^,F12.7,T,F6.2,'%n WRrTE(S,T INLET PROTEIN FLUX, W . Q C I n WRITE(S.T INLET PROTEIN FLUX, kgto * .QCin END IF IF (BCPSdn.Eai.DO) THEN WRITER-OUTLET PROTEIN FLUX, koto •'.QCout WRITE(B,TOUTLET PROTEIN FLUX, W . Q C o u t WRITER-OUTLET PROTEIN CONCENTRATION. kgftnS^.QCouVOoutS WRITE(S,*VOUTLET PROTEIN CONCENTRATION. kaftnS'.QCouVaoutS IF(HARV.E0.1)THEN WR(TE(S,TPROTEIN CONCN IN THE HARVESTING) RESERVOIR, k<ym3\ > (CTOTOCTOTI'VECSKrOTFLW WRITE(6,*)'PROTEIN CONCN IN THE HARVESTING RESERVOIR, kg/mS1, > (CTOT0-CTOT1*VEC8yTOTFLW END IF END IF WRITES/)" CTOT1,k0ft»S« ,.CTOT1 WRITE(S,T CTOT1, kotoiS - \CTOT1 IF (BCPSdn.Eai.D0.OR.BCPSup.Eai.D0) THEN WRITER,-)-AVERAGE CONCENTRATION, kg/mS *,CTOTAUVECS WRITE(S,*)'AVERAOE CONCENTRATION, kgVmS -'.CTOTAL/VECS END IF WRITE(8,7 AniTOT1,kg*.CTOT1 ,VEeS WRITE(6,T AmTOT1.kg*,CTOT1*VECS IF (BCP8dn.Eai .DaOR.BCPSup.Eai .D0) THEN WRITE(5,*)'TOTAL AMOUNT OF PROTEIN IN HFBR, kg *,CTOTAL WRITE(8,*)TOTAL AMOUNT OF PROTEIN IN HFBR, kg <CTOTAL IF(HARV.Eai )THEN WRIT^S.O'S PROTEIN REMOVED:',1.D2*(1.DO«TOT1*VECS/CTOT0) WRITES, Y H PROTEIN REMOVED:\1.D211.00-CTOTfVECS/CTOT0> END IF END IF C RETURN END C c.................~................................................. C SUBROUTINE SUBST(N,K,PS1,PS,C1,C) IMPLICIT REAL*a(A-H,L,aZ) DIMENSION PS1((fcN*1,0;K+1).P8«):NH,0:K*1) DIMENSION C1(0:N*1,0:K>1),C(0:N*1,0:K*1) C DO20M>,N+1 DO10J=O.K-H C1<I.J)-C(I,J) PS1(U)^S( I .J ) 10 CONTINUE 20 CONTINUE C RETURN END C c............................................ ..«.......««.. C SUBROUTINE OUTPUT <PS,PL.U,V,UL,CiX.RixOISP,RDI9P,CKAX,ITER,IS, > T,OT,NDISP,NDISPC,KDSP,KDISPC,INDC,INDPS,INOPo<,INDV) IMPLICIT REAL*8<A4l,L,OsZ> INTEGER EXF.ERF.ELF PARAMETER (KIS .N ' IOO) DIMENSION PL(0:N+1,0:K+1),PS(0:N-H,0:K*1).POlDSP(0:N,0:K) DIMENSION PLDSP(0:N.0:K),PSDSP(0-.N,0:K),PTOT(0:N,0:K) 0IMENSIONC(0:N*1,0:K+1),C0SP(0:N,0:K),CNSTRS(1:K) DIMENSION U(0:N,0:K+1).V(0:N*1.0:K),UL(0:N.0:K*1) DIMENSION UDSP(0:N,0:K),VDSP(0:N,0:K),ULDSP(0:N,0:K) DIMENSION XOISP(0:N),ROISP(0:K),X(0:N+1),R(0:K+1) COMMON /FI/FI.E.EL COMMON/PRESIT/IT COMMON /EXFRL/ EXF.ERF.ELF EXTERNAL FUNCTION FP C EX-OBLE(EXFyE ER=OBLE(ERFyE EXL4BLE(ELFyEL WRITE(6,S) WRITE(6,")TrERATION: '.ITER W R r T E f t T T . ^ . T . ' IS^.IS WRITE(6,*)'0T.»*,DT WRITE(S,*)'PRESSURE ITERATIONS: '.IT WRITE! S.YCMAX- ',CMAX WRITE(6.")-CMAX/0T- '.CMAX/DT S FORMAT^ C WRITER*) WRITE(6,408) <R0ISP<J)1.03,J=0,K,KDISPC) 408 FORMATf R,mni:'.F7.2,102F8.2) 40t FORMATfR.~™: ,,F8.2.102F«2) 410 FORMAT(-R.mm:',Fa.2,102F10.2) WRITES,*)' X,on CONCENTRATION FIELD (aohirf, intarpolMad):- DO27l=0,N*1 0O28JKP,K+1 IF (C(I.J).LT.0.DO) C(l,J)-0.D0 28 CONTINUE 27 CONTINUE CD8P(0,0)-C(0,0) CDSP(0,K>C(0,K+1) D O 1 0 I 1 . N - 1 CDSP(l,0H).6DO*(C(l'f1.0HC(l,0)) C0SP(l,K>0.6O0*(C(l+1,K>1)4C(l,K+1)) 10 CONTINUE CD3P(N,0)=C(N*1.0) CDSP(N,k>C(N+1,K>1) DO 26J*1.K-1 CD8P(0,J)=0.SOO*<G(0,J*1HC(0,J)) DO20l*1,N-1 CD3P(I.JH).2SO(r(C(I.JKC(l+1,J)*C(l,J*1>*C(l*1.J*1» 20 CONTINUE CD8P(N,JH).6O0-(C<N<-1.Jt1HC<N*1,J)) 26 CONTINUE DOS0NO.N.NDISPC WRITE(6,36)XDISP(l)*1.D2,(CDSP(l,J),J°0,K.KOISPC) SO CONTINUE 36 FORMAT(1X,F6.2,102F8.4) C IF( INDC.Eai)THEN WRITE(6,*) WR[TE(6.408) R<0).(R<J)1.D3,J=1.K*1,KDISPC) WRITE(6,<rX,eni CONCENTRATION FIELD aohul.aaiGulatad:' DO40l=O,N*1,NDI8PC WRITE<S,SS)X(I)1.D2,C(I,0),(G(I,J),J=1.K+1,KDISPC) 40 CONTINUE END IF C IF ( INDPOtEai .OR. INDP9.Ea i ) THEN D O 8 0 M . N DOSOJ-O.K PO«DSP(l,J)>*P(CDSP(M)) 80 CONTINUE • 0 CONTINUE END IF C IF(INDPS.Eai)THEN P8DSP(0,0)=PS<0,0) p808P(0,K)=P8(0,Kt1) 0O60MH.N-1 PSDSP(I.O)«0.5DO*<PS<I*1.0>4PS<I,0)> PSDSP<I,KH>.SD0"<PS<I*1,K*1>+PS<I,K*1)) 60 CONTINUE P8D8P(N.0)-P8(Nt1,0) PSDSP(N.K)H«(Nt1.K*1) D O 7 0 J 1 . K - 1 P3DSP(0,JH>.5DO*(PS<0.Jt1>+PS<0,J» D O 6 0 M . N - 1 PSD8P(U)=0.2600' > (PS(WHPS(WJH<>S(W+1HPS(I*1^+1)) 80 CONTINUE PS08P(N^)=0.600*(P8<N+1J+1)tP8(N*1.J)) 70 CONTINUE C D O S 1 M . N DO«1.M>,K PTOT<l,J)=PSDSP<l,JH>O»08P<l,J> 81 CONTINUE >1 CONTINUE C W R I T E R WRTTE(6,410) (RDI8P(J)*1.D V " * , K , K D 3 P ) WRITE(6,arX,oni ECS TOTAL PRESSURE FIELD:' OO100N0,N,NDISP WRITE<6,106)XDISP(I)1.02.(PTOT(UW=O.K,KDSP) 100 CONTINUE 106 FORMAT(F6.2,102F10.2) C W R I T E R WRITE(6.410) (R0I8P(JH.DS,J=0,K,KDSP) W R I T E R * . " * ECS HYDROSTATIC PRESSURE FIELD:' DO110l>O,N,NDISP WRITE(6,116)XDI8P(I)*1.D2,(P8D8P(I.J).JM),K.K08P) 110 CONTINUE 116 FORMAT(FS.2,102F10.2) END IF C IF(INDPOs.Eal)THEN WRITER W R I T E ( S , 4 1 0 ) ( R D I S P ( J ) - 1 . D 3 , J = 0 , K , K D 3 P ) W R I T E R ' X . E M E C S O S M O T I C P R E S S U R E F I E L D : ' DO117hO.N,N0ISP WRITE<6,1ie)XOI3P<iri.D2,(PO»D3P<l,J>,J=0,K.KDSP) 117 CONTINUE 118 FORMAT(F6.2.102F10.2) END IF C IF(INDPS.Eai)THEN PLD3P(0,0)=PL(0,0) PLD8P(0.KH>L(0,K*1) DO120l=1,N-1 PLD8P(I,0)=0.60(HPL(I+1.0)*PL(I,0)) PLDSP(I.KH>.S00"(PL(I-H,K+1)4PL(I,KH)) 120 CONTINUE PLDSP(N,0H>L(N«1,0) PLD8P(N,K)=PL(Nt1,K*1) DO-M0J-1.K-1 PLDSP<0,J)=O.6O0*<PL<0,J+1)*PL(0,J)) DO1S0M.N-1 PLD8P(l,J)=fl.2S0Cr(PL(UHPL(IH^HPL(I.Jt1)+PL(lt1.J*1)) 130 CONTINUE PLDSP(N.JH>.5O0*<PL(N*1.J+1HPL<N+1,J)) 140 CONTINUE C W R I T E R WRITE(6,410) (RDISP(J)1.DS,J=0,K,KD3P) W R I T E R X o m LUMEN PRESSURE FIELD:' DOISOMLN.NDISP WRrTE(5,15S)XDISP(l)1.D2,(PLD8P(l,J),J=0,K,KDSP) 160 CONTINUE 166 FORMAT(F6.2,102F10.2) END IF C UDSP(0.0)=U(0,0) ULDSP(0,OHJL(0,0) VO3P(0,0)=V(0,0) UDSP(N,OHI(N,0) ULDSP(N,OHIL(N,0) VD3P(N,0)=V(N*1.0) UDSP(0,K)=U(0,K+1) ULDSP(0,K>UL(0,K*1) VOSP(0.K)=V(0.K) UD3P(N,K)=U(N,Kt1) ULDSP(N,K)=UL(N,K»1) Appendix D: Source code in Fortran VD8P(N,K)*V(N<f1.K) DO160l*1,N-1 UDSP(l.a)»U(l,0) ULDSP(l,O)*UL(l.0) UDSP(I,K)*U(I,K*1) ULD8P(I,K>-UL<I,K*1> VDSF1I.OHI.DO V08P(l,KH).6O(nV(l,K)+V(l*1,K)) DO 105 J-1.K-1 UDSP<I,J)>0.600*(U(U)HI(M+1)) ULDSP(I,J)*).5D0>(UL(I,J)4UL<M+1)) VDSP<t,JH>.SO0-(V(l,J>tV<l*1,J)) 1SS CONTINUE 180 CONTINUE DO 170J1 .K -1 UDSP(O,J>*> .SOO* (U<O,JHU<O,J+ I ) ) U L D 8 P < 0 , J H > . S O 0 - < U M 0 , J > + U L < 0 , J + 1 ) ) VD8P<0.J)»V(0,J) UD3P(N,J)»O.SO01U<N,JHU<N,Jt1)> UL0SP(N,JH>.600*(UL<N,J>WJL<N.J+1)) VD8P(N,J)=V(N+1,J) 170 CONTINUE C W R ^ S . * ) WRrTE(6,40»MRDISP(J)*1.DWM),K,KD8P) WRITEfS,*) >'X,ora LUMEN AXIAL VELOCITY MhuJ,Mi rpoMad: UL'.ELF DO220M,N,NDISP WRITE(S.22S) XDISP(I)*1.D2.(UL0SPP^)«EXUJ>0,K,KDSP) 220 CONTINUE 226 FORMAT(FS.2,102F*.») C W R ^ S , * ) WRITE<5,40»XROISP(Jri.Daj*>,K,KDSP) WRITE(6,*) > ' X , c n ECS AXIAL VELOCITY »duil,lntapol»to<l: U-.EXF DO230l*>,N,N0ISP WRITE<5,235)X0ISP<l)1.D2.<U0SP<l,J)*ex,J=0,K.KDSP) 280 CONTINUE 23S FORMAT(F5.2,102F0.3) C WRrTE<S.*) WRITE(S,40tKR0ISP(J)*1.D3,J-4).K,K0SP) WRITE(5.T >'X.cra ECS RADIAL VELOCITY aokuMntirpoliUd: V . E R F DO240l»0.N,NDISP WRITE(5,24S)XDISP(I)*1.D2.(V0SP(U)*ER^'«.K,KDSP) 240 CONTINUE 246 FORMAT(F6.2,102F».4> C IF(INOV.EQ.1)THEN WRrTE<5,*) WRITE(5,400) R(0MR<J)*1.D3^*1,K*1,KDSP) WRITECS,") > 'X.on. LUMEN AXIAL VELOCITY aokid.edouMad: UL'.ELF 0O320l*>.N,NDISP WRITE(5,22S) > XDISP<I)"1.D2,UL<I,0)*EXL,<UL<I.J>*EXUJ=1,K+1,KDSP> S20 CONTINUE C WRITHS.") WRITE(6,4M) R<0).(R<J)*1.D3,J=1.K*1,KD3P) WRITE! 5,*) > 'X ,cn ECS AXIAL VELOCITY ackul.caleulatKl: U'.EXF 0O330M),N,NDISP WRITES, 235) XDISP(I)'1.D2,U(I,0)«EX,(U(I,J)*EX,J=1,K*1,KD3P) 330 CONTINUE C WRITER,") WRITE(6,40«) <RDISP<J)1.D3,J=0,K,KD8P) WRITHS,*) > •X.cni ECSRADIALVELOCrrYMihjd.ateutaUdiV.ERF DO 340 M.N+I .NOISP WRITE(5,245) X(P)1.D2.(V(I,J)*ER.J=0,K,KDSP) 340 CONTINUE END IF C RETURN END C c .................................................................... C SUBROUTINE INPF <IS.T.0T,ACCF,CTOTAL,CTOT0,TOTFLW.PL,PS,C) IMPUCIT REAL*8<A-H,L,»Z) PARAMETER (K«18, N=100) DIMENSION PL<0:N+1.0:K+1),PS(0:N<f1,0:K+1),C<0:N*1,0:K*1) C READ(UNrr"8)IS,T,DT,ACCF,CTOTAL,CTOT0,TOTFLW DO10l»0,N*1 READ(UNIT«3) (PL(I,J),J=0,K+1) 10 CONTINUE DO20M>,N*1 READ<UNIT*3) (PS<I,J),J=0,K+1) 20 CONTINUE DO30l=0,N+1 READ(UNITs3) (C(l J),J=O.K*1) 30 CONTINUE C RETURN END C c ...................................................... C SUBROUTINE OUTF <I8,T,DT,ACCF,CTOTAL,CTOTO,TOTFLW,PL,PS,C) IMPUCIT REAL*8tA-H,L,0-Z) PARAMETER (K*18,N=100) DIMENSION PL(O:N+1,0:K+1),PS<0:N+1,0:K-H>,C(0:NH,0:Kt1) C WRITER! NIT=4)IS,T.DT,ACCF,CTOTAL,CTOT0.TOTFLW D 0 1 0 l=O.N+1 WR[TE(UN(T=4) (PL(I.J).J=0,K*1) 10 CONTINUE DO20I-O.N+1 WRITE(UNfT=4) <PS<I,J),J=0,K*1) 20 CONTINUE DO30l=O,N*1 WRITE(UNIT=4) (C(l J),J=C,K*1) 30 CONTINUE B RETURN EN0 SUBROUTINE EXECTIME(HR1,HR2,MIN1,MIN2,8EC1,8EC2,HSEC1,HSEC2, > MOWTH1,MONTH2.DAY1,DAY2) INTEGER HR1,HR2,MIN1,MIN2,SEC1,SEC2,HSEC1,HSEC2 INTEGER MONTH1,MONTH2,DAY1,0AY2 INTEGER DHR,DMIN,DSEC.DHSEC,DMONTH,DDAY ' 0H8EC=HB6C2-HSEC1 IF(DHSEC.LT.0)THEN DHSEC*0H8EC+100 8EC2=8EC2-1 END IF D8EC=SEC2-SEC1 IF(DSEC.LT.0)THEN DSEC4SEC4O0 MIN2=*HN21 END IF DMIN4HN24MN1 IF(DMIN.LT.0)THEN DM1N=OMIN«0 HR2=HR2-1 END IF DHR°HR2-HR1 DMONTH4K>NTH24K>NTH1 DDAY=OAY2-0AY1 IF (DMONTH.Eai) THEN IF(MONTH1.Eai.OR.MONTH1.EaS.OR.MONTH1.Ea6.0R.MONTH1.EQ.7. » OR.MONTH1.Ea8.OR.MONTH1.EaiO.OR.MONTH1.EQ.12) THEN 0DAY=OAY2+31-OAY1 ELSE IF (WONTH1.EQ.2) THEN DDAY4AY2+284AY1 ELSE DDAY=OAY2+3<W>AY1 END IF END IF DHR=-OHR40DAY*24 WRfTE<6,10) DHR,DMIN,DSEC,DHSEC WRITE(9,10> DHR,DMIN,DSEC,DHSEC 10 FORMAT(/1X,'EXECUTION TIME: MS," h',13,' min',13,' •M4,' hs1) RETURN END
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Ukraine | 14 | 0 |
China | 10 | 1 |
United States | 5 | 0 |
Japan | 5 | 0 |
Russia | 3 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 16 | 0 |
Beijing | 10 | 0 |
Tokyo | 5 | 0 |
Tver | 3 | 0 |
Redmond | 1 | 0 |
Mountain View | 1 | 0 |
Ashburn | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: