Computational and Mathematical Methods in Medicine

Volume 2012 (2012), Article ID 973424, 16 pages

doi:10.1155/2012/973424

Research Article

Finite Volume Scheme for Double Convection-Diffusion Exchange of Solutes in Bicarbonate High-Flux Hollow-Fiber Dialyzer Therapy

Kodwo Annan

Department of Mathematics and Computer Science, Minot State University, Minot, ND 58707, USA

Received 20 May 2012; Accepted 8 August 2012

Academic Editor: Timothy David

Copyright © 2012 Kodwo Annan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The efficiency of a high-flux dialyzer in terms of buffering and toxic solute removal largely depends on the ability to use convection-diffusion mechanism inside the membrane. A two-dimensional transient convection-diffusion model coupled with acid-base correction term was developed. A finite volume technique was used to discretize the model and to numerically simulate it using MATLAB software tool. We observed that small solute concentration gradients peaked and were large enough to activate solute diffusion process in the membrane. While CO₂ concentration gradients diminished from their maxima and shifted toward the end of the membrane, concentration gradients peaked at the same position. Also, CO₂ concentration decreased rapidly within the first 47 minutes while optimal concentration was achieved within 30 minutes of the therapy. Abnormally high diffusion fluxes were observed near the blood-membrane interface that increased diffusion driving force and enhanced the overall diffusive process. While convective flux dominated total flux during the dialysis session, there was a continuous interference between convection and diffusion fluxes that call for the need to seek minimal interference between these two mechanisms. This is critical for the effective design and operation of high-flux dialyzers.

1. Introduction

High-flux dialyzer is one of the possible treatments to remove toxic solutes from the blood when the native kidneys loss their function. Small solutes removal is primarily done by diffusion while larger solutes removal is obtained by convection. The efficiency of a dialyzer is therefore dependent on its ability to use these mechanisms (convection and diffusion) to exchange solutes across the dialyzer membrane [1–4]. Diffusion is mainly affected by blood and dialysate flow rates, dialyzer surface area, temperature, and membrane thickness. If we assume constant values to all other factors, then the diffusion mechanism depends on the blood and dialysate concentration gradients [5, 6]. This, however, is influenced by the blood and dialysate flow distributions and flow rates. Extensive research has been done on flow distribution mismatch frequently observed at the blood-dialysate interface [6–10]. Physically, attempts to correct and optimize blood and dialysate flow mismatch have been made by redesigning blood and dialysate headers. Options such as space yarns and moiré structure have been proposed to resolve dialysate channeling phenomenon external to the fiber bundle [11, 12]. The main feature of convection is the use of high-flux HD characterized by high permeability for water, electrolytes, and higher clearance of middle and large molecular weight solutes. The role of convective transport is discussed extensively in recent articles [13–19].

The investigation of the effect of convection and diffusion during dialysis session continues to pose a major challenge to HDF researchers and engineers. Best-known earlier models were based on clinical data. However, these macroscopic experimental approaches make it difficult to capture and explore convective and diffusive transports during dialysis session. Mathematical models have been used to evaluate, optimize, and control various forms of dialysis therapy from clinical routine to investigating new issues in dialysis therapy [3–6, 8–11, 15, 19, 20]. The underlying mechanism of these mathematical models has been Navier-Stoke equation. Numerical methods aimed to quantify both the convective and diffusive transports of solutes exchange across membranes have been used [21–28]. Researchers [21–23, 25, 27, 28] used finite difference schemes and control volumes while analytical solutions were derived by [24, 26]. These authors, however, neglected the effects of either diffusion or convection flows, or their choice of techniques not justified in dialysis therapies, or did not include buffer which is common in dialysis sessions.

In dialysis, the type of numerical scheme used in computing solutions of convective-diffusive equations is very necessary, especially when high flux membrane (i.e., where convection term dominates) is used. In this case, the dialyzer membrane is so thin that one is forced to use under resolved methods that may be unstable. On the other hand, the use of dispersive schemes may trigger numerical instabilities which may affect the fully description of the convection and diffusion phenomena during HD session. To achieve maximal dialyzer efficiency, the accuracy and reliability of numerical schemes used to compute convection-diffusion phenomena is of paramount importance. These numerical schemes depend on the choice of discretization and the quality of the underlying mesh.

This paper focused on finite volume method (FVM) for unsteady-state convection-diffusion equations that arise in dialysis therapy. The transport equation was described using a three-compartmental model of blood, membrane, and dialysate compartments. The model was coupled with buffering and replenishment. An accurate transient convection-diffusion model that described solute exchange in a typical high-flux hollow-fiber dialyzer was performed. The numerical discretization and analysis schemes were then proposed and tested. We then explored the impact of small molecule weight solute (carbon dioxide and bicarbonate) transports in a high-flux dialyzer followed by conclusions.

2. Model Formulation

2.1. Membrane Model

The conservation law for the transport of solute concentrations in an unsteady flow has the general form

where (a constant) is the density of incompressible fluid, is the solute concentration, is the fluid velocity, is the production of new solute at that point, is the diffusion constant, and and are the normal vector operators. Equation (1) basically states that the rate of increase of the number of molecules of solute () at any point equals the (negative of the) rate they are being removed at that point by convection () plus the rate they are being added by diffusion () plus the rate at which solutes are being produced ().

The following simplifying assumptions are made in the membrane model.

(1)The membrane is small enough that it is assumed to be in equilibrium (steady state), so the time derivative term is zero.(2)The membrane impedes flow in all directions but radially, so the velocity vector is in the direction only.(3)The membrane volume is small that production of new solute can be ignored, so .(4)Since we are interested in the change of concentration from one side of the membrane to the other, the rate of change of concentration in the axial direction is smaller (and it is zero in the direction due to symmetry). Therefore, we ignore the terms in except for the radial direction. Using assumptions 1–4, (1) becomes

Integrating (2) and transforming the resulting equation using the divergence theorem gives

where the normal vector to the face is , and the integrals indicate either a volume integral over the control volume (with ) or a surface integral over the faces (with ). All components of in (3) are zero except the radial direction, so the dot product with has only the term . Using assumption 4 and assuming that the functions in the integrands are constant across the faces that the surface integrals act on, (3) reduces to

From (4), if the fluid velocity was zero, the LHS would be zero, the gradient would be a constant, and the concentration would follow a linear slope through the membrane in the radial direction. However, with a nonzero velocity, the flux of solute due to convection, , is nonzero. Since the concentration varies through the membrane in the radial direction, intuitively, the diffusion term would need to compensate for the drop in convection so that a constant flow exists across the membrane from one side to the other. Re-arrange the terms in (4)

Since (5) is true for any separation and for any value , the total solute flux, (both convection and diffusion) is constant across the membrane at any point. Thus, we have

The exact solution for (6) in terms of the concentration is of the form

Substituting (7) into (6) and gathering terms

The coefficient of the exponential must be zero and the constant term equal , thus

The constant is determined by the concentration value at one end of the membrane. For this paper, the pressure in the dialysate side (at larger values of ) is normally larger than that in the blood side, so the velocity is negative. Picking a solute where the concentration is higher in the dialysate side, gives positive concentration gradient, . Therefore, if both terms of (6) are negative, then is negative, that is, a total flux in the negative direction toward the blood side. In this case, is positive, is negative, and for to be positive must be negative. So the concentration function must be of the form

where the ’s are the positive versions of the constants. Therefore, concentration is positive and decreasing for smaller (toward the blood side) and the slope is increasing for smaller to compensate for the reduced convection.

2.2. Transmembrane Flow

Following [29] and assuming that reflection coefficient is negligible because of small () fiber pore size [30], we describe the flow passing through the membrane (see Figure 1(b)) by simplified Kedem-Katchalsky (K-K) equations

(m/s) is ultrafiltration velocity or volumetric flux across the membrane; is solute flux across the membrane; is the hydraulic permeability of the membrane; is solute diffusive permeability coefficient of a membrane; represents the average solute concentration at each side of the membrane; is solute concentration difference (i.e., transmembrane concentration) across the membrane. The parameter is the membrane surface hydraulic permeability of the membrane. Thus, the membrane interfacial conditions for the blood-side model are

(a)

(b)

Figure 1: (a) Schematic of a typical hollow-fiber dialyzer module with the computational blood and dialysate domains. (b) Mass transport of solutes in blood and dialysate compartments through a single hollow fiber membrane.

2.3. Blood-Side Flow Model

Consider () as coordinates representing a point in the cylindrical coordinate system where the -axis is taken along the dialyzer length (i.e., ) and is taken along the radial direction. An axisymmetric domain, where is chosen to lie in the range between and for a membrane length and radius   (see Figure 1) depicts the blood side model. The Navier-Stokes and continuity equations that govern the flow of an incompressible Newtonian fluid representing blood with constant density and viscosity can be described as [13]

where and are the radial and axial velocity components, respectively, and the pressure. Using the continuity equation and the fact that flow is driven by pressure gradient in the -direction, a fully developed inlet velocity profile for number of fibers at and are obtained [31, 34]

Here, is the inlet blood flow rate in each of the hollow fibers with a fiber cross-section area . Applying no slip condition at the wall and axisymmetric axis, respectively, at

The convection-diffusion equation governing the mass transport of solutes coupled to the blood velocity field is given by

where and are the concentration and the diffusion coefficient of solute in the blood, respectively. The inlet and outlet boundary conditions for the concentration equation (4) are

defines the buffer term that vanishes everywhere except in the blood membrane domain and denotes the rate of solute production or consumption per time. We adapt buffer reaction rates for given by [31, 35]

where ,  ,   and  / = 20. The rate controlling reactions for the carbonate and bicarbonate ions are given as [31, 35]

where , and are their reaction constants with their equilibrium constants defined as and , respectively. The overall reaction is

with the following fast reactions assumed to be at equilibrium

where and are their equilibrium constants. The parameters and their values are stated in Table 1.

Table 1: Reaction and equilibrium constants and equations used in this paper at 297 K.

Constant
Value/equation
Unit
Ref

Forward reaction constant,
2.38 × 10ࢤ2

[31, 32]

Reverse reaction constant,
1.4 × 10¹

[31, 33]

Forward reaction constant,

[31, 32]

Reverse reaction constant,

[31, 32]

Equilibrium constant,

[32]

Equilibrium constant,

[32]

Equilibrium constant,

[32]

Equilibrium constant,

[32]

2.4. Dialysate-Side Flow

Since each fiber was surrounded by a uniform annulus (shown in Figure 2(a)), we adapted Krogh cylinder geometry [36, 37] with annulus radius which was far larger than the fiber radius . Assuming a fully developed axial and radial velocities in annulus geometry, the generic continuity and momentum equations reduced to (22) as reported by [38], with specified boundary conditions of at and

Here, the parameter represented flow rate in the dialysate inlet, the ratio of , and the width of the raised collar used to promote uniform flow in dialyzers. The solute replenishment term, , was introduced to help maintain dialysate concentration level and was calculated using [31]

where is the replenishment coefficient.

(a)

(b)

Figure 2: Geometric configuration near dialysate inlet. (a) The distributors are designed to keep the dialysate entering the hollow-fibers uniformly. (b) Schematic computational domain at the dialysate compartment.

The transport of solutes in the annulus, shown in Figure 2(b), involving convection and diffusion with   and   defined by (22) could be described similarly as (16) as

3. Algorithm and Numerical Techniques

Finite volume method (FVM) was used to transform the model equations (12)–(24) into dimensionless system. Since the structure of the convection-diffusion equations (16) and (24) only differed by the source term, we replaced the source term by .

3.1. Transformation of Models Using FVM

Integrating both sides of (16) or (24) over a small control volume gave

Thus, (25) means that the rate of increase of concentration with time in the volume element is equal to the convective flow into the volume element, plus the diffusive flow, and the creation of new solute from the source term totaled over the volume element. Since both convective and diffusive terms represented divergence of vector fields (i.e., the fluid flow vector and the concentration gradient, resp.), we applied the divergence theorem to the integrals of these terms to convert them to surface integrals.

3.1.1. Convective Term

Since the flow velocity is only in the -direction, the convective flow is

Therefore, the divergence of the convective flow vector using cylindrical coordinate is

where we have used the fact that the flow components in the and directions are zero and that the component is a function of only, implying . Thus, the second term on the left-hand side of (25) is

Since the vector is in the -direction, it does not cross the surfaces of the control volume cube on the faces in and directions. That is, the normal to those faces are perpendicular to the flow vector and so the dot product is zero. Therefore, the only nonzero parts of the surface integral are those over the faces in and directions of the cube. The normal to those faces is parallel to the convection vector (in the direction and antiparallel in the direction), so the integral becomes

where is the size of the volume cube in the -direction and is the coordinate at the center of the cubic volume. In the process of discretization, we approximate the values of and over the surface area of the cube by their values at the nearby grid points as and , respectively, if the grid point is sufficiently fine. Thus,

3.1.2. Diffusive Term

Since diffusion is driven by concentration gradient, the divergence vector field in cylindrical coordinates is

Using the divergence theorem, the diffusion term in (25) could be written as

For the surface integral on the and faces, only the first element of the gradient vector is applicable (the normal to those faces picks out that component of the vector) while the second element is used on the and faces only. As a result, the integral (32) becomse

Approximating the values of and over the face area by indicating their values with subscript “” and pull the constant values out of the integral, right-hand side of (33) becomes

Thus, the LHS of (25) with defining the volume average of solute concentration is

The RHS of (25) becomes

where .

3.2. Scaling to Dimensionless Form

The transformed equations (35) and (36) and their initial and boundary conditions (14)-(15) and (17)–(23) were converted into nondimensional forms using the same scale factors for both blood and dialysate flow regions. The nondimensional variables used in the transformation are indicated with superscript “*” below and the reference variables defined in Table 2

where and are Pe’clet and Sherwood numbers, respectively, and the ratio of momentum diffusivity and mass diffusivity is donated by . expressed the relative importance of convection to diffusion while related inertial effects to viscous effects. Since dialysis devices employ laminar fluids flow with the inertial effects would be irrelevant [31, 41].

Table 2: The reference variables with their description.

Symbol
Description

Reference length: they are the same for both compartments

Reference velocity

Reference solute concentration

Substituting the dimensionless variables in (37) into (35) and (36), simplifying notations, and dropping the superscript “*” resulted in

The dimensionless initial and boundary conditions, buffer and replenishment, and membrane interfacial conditions corresponding to (38) were as follows.

Blood-Side Inlet Velocity Conditions

Dialysate-Side Inlet and Outlet Velocities

No Slip and Axisymmetric Conditions

Inlet and Outlet Blood and Dialysate Concentrations

Buffer and Replenishment Terms

where represented dimensionless buffer terms in blood side and depicted dimensionless replenishment term for solute and .

Blood-Membrane Interfacial Conditions

where is the Sherwood number and is the ratio of momentum and mass diffusivity defined in (37).

3.3. Model Parameters and Numerical Algorithm

3.3.1. Geometric and Transport Parameters

The hollow-fiber dialyzer chosen for this study was the Fresenius’ F60 model with membrane area was 1.15 m². The membrane module has 22 cm effective axial length with 200 μm and 40 μm fiber diameter and thickness, respectively [20]. The initial inlet bicarbonate concentration values of blood and dialysate were set to 19 mol·m⁻³ and 35 mol·m⁻³, respectively, while the blood-side and dialysate-side flow rates were, respectively, 400 mL/min (i.e.,      m³ s⁻¹) and 800 mL/min ( m³ s⁻¹) [31, 42]. Other parameters and constant values used in this paper are either listed in Table 3 or are computed using values in Table 3.

Table 3: Geometric and transport characteristics of the hollow-fiber module used.

Parameter (unit)
Notation
Value
Ref.

Diffusion coefficient of in blood ()

[31–33, 35]

Diffusion coefficient of in blood ()

[31–33, 35]

diffusion coefficient in dialysate ()

[31–33, 35]

diffusion coefficient in dialysate ()

[31–33, 35]

Membrane effective length (m)

0.20
[31]

Hydraulic permeability ()

[31]

Width of raise collar (m)

0.014
[31]

Fiber diameter (μm)

200
[31, 39]

Fiber thickness (μm)

40
[40]

Number of fibers

9000–12000
F60 Model

Membrane permeability of ()

[32, 33]

Membrane permeability of ()

[32, 33]

Radius of dialysate channel (m)

[39, 40]

Radius of blood channel (m)

[39, 40]

Initial velocity at blood inlet ()

[31, 39, 40]

Initial velocity at dialysate inlet ()

[31, 39, 40]

3.3.2. Variables and Grid Definition

Application of FVM resulted in the creation of grid structured such that the number of rectangular cells in and direction remained constant throughout the domain of interest. For the spatial domain, the numerical model used separate subdomain grids for the blood side and the dialysate side since the two models and their domain dimensions were different. The spatial grid had a variable number of intervals in each axis, defined by the variables , , and . Because of the FVM development, the boundaries of the domain of interest have to be at the faces of the rectangular control volumes, rather than at a grid point. Thus, for example, the inlet boundary condition applied at , so the first internal grid point is at , where defined the grid spacing. However, one extra grid row or column outside the boundary was used to allow easy application of boundary conditions. Thus, the first grid point in was outside the inlet at . Therefore, the grid size in each axis was

since the domain of interest in the model was for the -direction and and in the -direction for blood and dialysate sides, respectively.

The indices of the variables (such as solute concentration) ran from to and to or . Therefore, the boundaries of the spatial domains were between indices and 2 and between and (blood side), similarly for the dialysate side, it was between and 2 and between and . Since the boundary conditions define the values of the variables there, the numerical model only calculated the values in the range from 2 to and so on.

For the time coordinate, a time increment of was used to sample the time axis. The real time represented by a sample is . Finally, we considered two types of solutes defined by the subscript in the models, as per the Table 4.

Table 4: Solutes indexes for numerical computations.

Value for s index
Solute type

1
CO₂ (carbon dioxide)

2
(bicarbonate)

Therefore, the variables in the domain of interest could be represented by a 4-dimensional array where the variable can be , and .

3.3.3. Hybrid Differencing Scheme

The continuous diffusion terms in (38) were numerically discretized using the hybrid differencing scheme. The scheme switched to the upwind differencing when the central differencing produced inaccurate results at high Peclet numbers. Also, since the partial continuous derivatives were at the faces of the control volume, we used the values at the center and adjacent grid points to find the values of and . The convection terms in (38) used the upstream scheme to estimate the value of the function at the control volume face. Thus, the functional values at the faces of our model system were represented by the following.

Blood Side:

Dialysate Side:

3.3.4. Boundary Conditions and Stability

Since many of the boundary conditions (BCs) depended on values in adjacent grid points, and the numerical equations for our model system only defined the values in the interior of the domain, the BCs were defined in terms of grid values to the interior of that boundary. Therefore, the corners were resolved by simulating the interior points of the model system followed by the BCs of the blood, dialysate, and the membrane (see Figure 3).

Figure 3: Dependencies between the BCs with arrows, which in turn defined the sequence of application required by these dependencies. The BC at the upstream boundary in each side was a constant and all other BCs depended on values in the current time step [31].

4. Results and Discussions

The numerical solution presented in the previous sections allowed us to determine the concentration gradients, convective and diffusive fluxes, total flux, and concentrations of carbon dioxide and bicarbonate profiles for high-flux membrane. Diffusive solute transport across the membrane was predominantly driven by concentration gradients, whereas convection transport was determined by pressure gradients.

4.1. Carbon Dioxide and Bicarbonate Concentration Gradients

Concentration gradient has been the principal process for removing end-products of metabolism (urea, creatinine, uric acid) and for repletion of bicarbonate deficit of metabolic acidosis associated with end-stage renal disease patients. Figure 4 showed the results of the concentration gradient profiles for carbon dioxide and bicarbonate solutes in the membrane for different time periods. The membrane carbon dioxide concentration gradient profiles (see Figure 4(a)) increased, reaching maxima gradients inside the membrane, then appeared to diminish from their maxima along the dialyzer axial length. The maximum points do shift toward the end of the membrane length while the maximum magnitude occurred at minutes.

(a) Concentration gradients at various times

(b) Concentration gradients at various times

Figure 4: Bicarbonate and carbon dioxide concentration gradients increased along the membrane axial distance at different time periods. (a) Depicted carbon dioxide maxima gradients shifting toward the membrane end and decreasing from their maxima. (b) Displayed bicarbonate gradients peaking at the same position inside the membrane. Maxima concentration gradients were achieved for both carbon dioxide and bicarbonate under an hour during the dialysis therapy.

Figure 4(b) depicted positive increased bicarbonate concentration gradients at different time periods in the membrane as the axial distance increased. The gradients peaked around the same position at and appeared to experience reduction of the bicarbonate concentration gradients toward the end of the membrane. Since bicarbonate containing dialysate was used in our model, it was important to have adequate concentration gradients to generate bicarbonate flux into the blood to restore body buffer. The adequacy of bicarbonate concentration gradient was observed as shown in Figure 4(b).

4.2. Carbon Dioxide and Bicarbonate Diffusive Fluxes

The fluxe profiles for carbon dioxide and bicarbonate in the membrane were presented at the membrane region ( cm) for different time periods in Figure 5. Both 5(a) and 5(b) showed the unsteady characteristics of solute diffusive fluxes at various radial positions in the membrane. In all the chosen radial locations, the rate of mass transfer of the and solutes increased at the onset of the dialysis therapy followed by a small fluctuations and then became constant during the remaining therapy session. In Figure 5(a), the sharp increased in may be caused by (i) the active hydrogen ion from the blood reacting with bicarbonate ions to produce more of the and/or (ii) the incomplete dissociation of bicarbonate ion into in the membrane. Similarly, the initial bicarbonate sharp increase observed may be explained by the incomplete dissociation of carbonic acid into bicarbonate and hydrogen ions. These observations suggested a bicarbonate ion carryover effect in the membrane during dialysis therapies. In addition, both figures depicted increased diffusive fluxes toward the membrane end and an abnormally high solute fluxes near the blood-membrane interface. This abnormality may increase the driving force for diffusion and eventually enhance the diffusive process during dialysis session.

(a) Diffusion fluxes at various radial positions

(b) Diffusive fluxes at various radial positions

Figure 5: Bicarbonate and carbon dioxide diffusive fluxes in the membrane at dialyzer axial distance  cm. Diffusive solute fluxes increased as the radial distance increased. Also, diffusion for both solutes was stronger near the fiber walls than the center. The carbon dioxide and bicarbonate ions effect observed before solute stability could enhance the diffusion process during dialysis session.

4.3. Carbon Dioxide and Bicarbonate Convective Fluxes

Figure 6 showed the results of carbon dioxide and bicarbonate convective solute fluxes in the dialyzer membrane for various radial positions in the first hour of dialysis session. Compared to Figure 5, Figure 6 displayed the dominance of convective fluxes for small solutes ( and ) during a high-flux hollow-fiber dialysis session. The magnitude of and fluxes were determined by the natural convection (transmembrane pressure gradient) and the forced convection (mass inflow). Both figures also indicated that convective fluxes decreased in the membrane with increased radial distance as one moved toward the end of the axial length. Thus, in this model there was a continuous interference between convective and diffusive fluxes (see Figures 5 and 6). In this case, increasing one type of transport mechanism would decrease the other and therefore could be beneficial or detrimental on dialyzer’s efficiency. In the region near the blood ports, the sagging nature of the convective profiles may explain the abnormally high diffusive flux of solutes at the blood-membrane interface observed in Figure 5. This observation craves the need to seek minimal interference between convection and diffusion during dialysis therapy.

(a) convective fluxes at various radial positions

(b) convective fluxes at various radial positions

Figure 6: Various carbon dioxide and bicarbonate convective fluxes at different radial distances were presented. Contrary to diffusive flux, convective flux decreased at the fiber wall and increased in the fiber center. Also, the sagging nature of convective fluxes observed at the blood port for the solutes may explain the high diffusive flux abnormalities observed in Figure 5.

4.4. Total Fluxes in Dialysis Membrane

Figure 7 displayed the dominance of convective flux over diffusion when a high-flux hollow-fiber dialyzer was used. Both figures showed that the total flux (convection and diffusion) of solutes were mediated by convective flux and it decreased along the axial length. The decreasing profile within the membrane may be explained by the decreasing nature of the transmembrane pressure gradient or the solute accumulation at the membrane surface over time. Thus, in addition to convection playing a major role in higher molecular weight solute transports [13], it was also more efficient than diffusion in small solutes transport when high-flux dialysis membrane was used.

(a) Convective and total fluxes in membrane

(b) Convective and total fluxes in membrane

Figure 7: Total fluxes for both carbon dioxide and bicarbonate solutes were mediated by convective fluxes. Total fluxes decreased along the dialyzer distance.

4.5. Carbon Dioxide Concentration

Carbon dioxide concentration profiles in the membrane at various time periods and radial distances during high-flux dialysis session were shown in Figure 8. In Figure 8(a), various concentration profiles at different time periods were shown as a function of membrane axial distance. The result clearly showed that the concentration decreased as the axial length increased and that most desorption occurred within 47 minutes of the therapy session. The constant profiles after 47 minutes indicated that the membrane was fully saturated with concentration during the therapy. In Figure 8(b), concentration profiles as a function of the first 70 minutes of dialysis session was shown at different radial positions. The result demonstrated that concentration increased as one moved toward the end of the membrane. Also, the concentration near the wall of the fiber was much higher than that of the center of the fiber at the same axial position. At the outlet, the concentration was stable with respect to radial position, indicating saturation in the membrane.

(a) Concentrations at different time periods

(b) Concentrations different radial distances

Figure 8: The effects of diffusion and convection fluxes and concentration gradients for carbon dioxide solute at different time periods and radial distances were presented. Carbon dioxide decreased over time implying carbon dioxide solute desorption in the membrane (a). Also, carbon dioxide concentration was higher in the fiber center than near the fiber walls (b).

4.6. Bicarbonate Concentration

Similarly, bicarbonate concentration profiles in the membrane at various time periods and radial positions during high-flux dialysis session were shown in Figure 9. At various time periods, the result indicated a concentration increased as the membrane axial distance increased (see Figure 9(a)). Also, it was shown that concentration increased rapidly within the first 30 minutes and then became stable after 40 minutes of the therapy session. In Figure 9(b), the concentration profiles at axial distance (m) against time for various radial distances was shown. concentration increased and was nearly saturated in the membrane. The rate of increased and the degree at which concentration increased may be determined by the immediate buffer response through convection and diffusion and the extent to which organic acid production in the membrane is increased. In addition, these observations confirmed clinical studies [31, 42] that dialysis patients achieve stable physiologic concentration levels during dialysis therapy.

(a) Concentrations at various time periods

(b) Concentration for various radial positions

Figure 9: Bicarbonate concentration profiles at various time periods and radial distances in the membrane were presented. At different time periods, bicarbonate concentration increased rapidly for during the first 30 minutes and then slowly for another 10 minutes, before becoming stable for the rest of the therapy (a). The maximum concentration reached was within the physiologic range reported by clinical studies. In (b), the bicarbonate concentration profiles within the first 40 minutes at various radial distances were found to increase until the membrane becoming nearly saturated with the solute. However, concentration at the fiber walls was much higher than that in the fiber center.

5. Conclusions

A mathematical model that coupled nonlinear unsteady convection-diffusion mass transfer of small solutes in high-flux dialyzer with buffer during dialysis session was developed. Finite volume technique was used to transform the model equations and numerical discretization and analysis schemes were then proposed and tested. The solute concentration gradients, diffusive and convective fluxes, and their effects on overall concentrations were explored. The salient observations were summarized as follows:

 Both carbon dioxide and bicarbonate concentration gradients increased and reached maxima gradient values inside the dialyzer membrane. While concentration gradients appeared to diminish from their maxima and shift toward the end of the membrane, concentration gradients peaked at the same position. In addition, the magnitude of the concentration gradient was large enough to activate diffusion in the membrane.  Diffusive fluxes for carbon dioxide and bicarbonate showed increased profiles at different radial distances in the membrane. Abnormally high fluxes were observed near the blood-membrane interface that could increase diffusion driving force and eventually enhance the overall diffusive process during dialysis session. Convective flux still dominated total flux for small solute ( and ) transfer during high-flux dialysis therapy. However, the continuous interference between convective and diffusive fluxes could be beneficial or detrimental when accessing high-flux dialyzer efficiency for small solute transport. Carbon dioxide concentration decreased rapidly causing the membrane to become fully saturated with within 47 minutes. Further investigation showed an increased concentration towards the end of the membrane which indicated higher concentration near the walls of the fiber than the fiber center at the same axial distance. Similarly, there was an optimal concentration in the membrane within 30 minutes of dialysis therapy because of the effects of convection and the preponderance of diffusion. The rate of increased and the degree at which bicarbonate increased could be caused by immediate buffer response and the positive effect of convection on diffusion or the extent to which organic acid production in the membrane was increased.

Therefore, the model presented provided an accurate quantitative description of both convection and diffusion through high-flux membrane with buffer. This is critical for the effective design and efficient operation of dialyzers.

References

K. Sakai, “Technical determination of optimal dimension of hollow fibre membranes for clinical dialysis,” Nephrology Dialysis Transplantation, vol. 4, supplement 6, pp. S73–S77, 1989.

J. E. Sigdell, “Calculation of combined diffusive and convective mass transfer,” International Journal of Artificial Organs, vol. 5, no. 6, pp. 361–372, 1982.

C. Ronco, A. Brendolan, C. Crepaldi, M. Rodighiero, and M. Scabardi, “Blood and dialysate flow distributions in hollow-fiber hemodialyzers analyzed by computerized helical scanning technique,” Journal of the American Society of Nephrology, vol. 13, no. 1, pp. S53–S61, 2001.

C. Ronco, B. Breuer, and S. K. Bowry, “Hemodialysis membranes for high-volume hemodialytic therapies: the application of nanotechnology,” Hemodialysis International, vol. 10, no. 1, pp. S48–S50, 2006.

C. V. Velde and E. F. Leonard, “Theoretical assessment of the effect of flow maldistributions on the mass transfer efficiency of artificial organs,” Medical and Biological Engineering and Computing, vol. 23, no. 3, pp. 224–229, 1985.

N. Hoenich and C. Ronco, “Dialyzers evaluation,” in Replacement of Renal Function by Dialysis, J. Winchester, K. Koch, C. Jacobs, and C. Kjiellestrand, Eds., pp. 256–270, Kluwer Academic, Dordrecht, The Netherlands, 1996.

S. Takesawa, M. Terasawa, M. Sakagami, T. Kobayashi, H. Hidai, and K. Sakai, “Nondestructive evaluation by X-ray computed tomography of dialysate flow patterns in capillary dialyzers,” ASAIO Transactions, vol. 34, no. 3, pp. 794–799, 1988.

C. V. Velde and E. F. Leonard, “Theoretical assessment of the effect of flow maldistributions on the mass transfer efficiency of artificial organs,” Medical and Biological Engineering and Computing, vol. 23, no. 3, pp. 224–229, 1985.

K. Annan, “Mathematical modeling for hollow fiber dialyzer: Blood and dialysate flow characteristics,” International Journal of Pure and Applied Mathematics, vol. 79, no. 3, pp. 425–452, 2012.

A. L. Babb, P. C. Farrell, D. A. Uvelli, and B. H. Scribner, “Hemodialyzer evaluation by examination of solute molecular spectra.,” Transactions of the American Society for Artificial Internal Organs, vol. 18, no. 0, pp. 98–122, 1972.

F. Gastaldon, A. Brendolan, C. Crepaldi et al., “Effects of novel manufacturing technology on blood and dialysate flow distribution in a new low flux “α polysulfone” hemodialyzer,” International Journal of Artificial Organs, vol. 26, no. 2, pp. 105–112, 2003.

C. Ronco, S. Kumar Bowry, A. Brendolan et al., “Hemodialyzer: from macro-design to membrane nanostructure; the case of the FX-class of hemodialyzers,” Kidney International, vol. 61, supplement 80, pp. S126–S142, 2002.

F. Locatelli, C. Manzoni, and S. Di Filippo, “The importance of convective transport,” Kidney International, vol. 61, supplement 80, pp. S115–S120, 2002.

M. Y. Jaffrin, L. Ding, and J. M. Laurent, “Simultaneous convective and diffusive mass transfers in a hemodialyser,” Journal of Biomechanical Engineering, vol. 112, no. 2, pp. 212–219, 1990.

J. E. Sigdell, “Calculation of combined diffusive and convective mass transfer,” International Journal of Artificial Organs, vol. 5, no. 6, pp. 361–372, 1982.

M. Y. Jaffrin, “Convective mass transfer in hemodialysis,” Artificial Organs, vol. 19, no. 11, pp. 1162–1171, 1995.

A. Werynski and J. Waniewski, “Theoretical description of mass transport in medical membrane devices,” Artificial Organs, vol. 19, no. 5, pp. 420–427, 1995.

A. Wüpper, D. Woermann, F. Dellanna, and C. A. Baldamus, “Retrofiltration rates in high-flux hollow fiber hemodialyzers: analysis of clinical data,” Journal of Membrane Science, vol. 121, no. 1, pp. 109–116, 1996.

J. M. Miranda and J. B. L. M. Campos, “An improved numerical scheme to study mass transfer over a separation membrane,” Journal of Membrane Science, vol. 188, no. 1, pp. 49–59, 2001.

C. Legallais, G. Catapano, B. V. Harten, and U. Baurmeister, “A theoretical model to predict the in vitro performance of hemodiafilters,” Journal of Membrane Science, vol. 168, no. 1-2, pp. 3–15, 2000.

V. Calabrò, S. Curcio, and G. Iorio, “A theoretical analysis of transport phenomena in a hollow fiber membrane bioreactor with immobilized biocatalyst,” Journal of Membrane Science, vol. 206, no. 1-2, pp. 217–241, 2002.

V. Calabrò, S. Curcio, and G. Iorio, “A theoretical and experimental analysis of a membrane bioreactor performance in recycle configuration,” Journal of Membrane Science, vol. 273, no. 1-2, pp. 129–142, 2006.

N. S. Abdullah, D. R. Jones, and D. B. Das, “Nutrient transport in bioreactors for bone tissue growth: why do hollow fibre membrane bioreactors work?” Chemical Engineering Science, vol. 64, no. 1, pp. 109–125, 2009.

C. Heath and G. Belfort, “Immobilization of suspended mammalian cells: analysis of hollow fiber and microcapsule bioreactors,” Advances in Biochemical Engineering/Biotechnology, vol. 34, pp. 1–31, 1987.

M. Labecki, B. D. Bowen, and J. M. Piret, “Two-dimensional analysis of protein transport in the extracapillary space of hollow-fibre bioreactors,” Chemical Engineering Science, vol. 51, no. 17, pp. 4197–4213, 1996.

K. Li and X. Tan, “Mass transfer and chemical reaction in hollow-fiber membrane reactors,” AIChE Journal, vol. 47, no. 2, pp. 427–435, 2001.

A. Y. Patkar, J. Koska, D. G. Taylor, B. D. Bowen, and J. M. Piret, “Protein transport in ultrafiltration hollow-fiber bioreactors,” AIChE Journal, vol. 41, no. 2, pp. 415–425, 1995.

P. M. Salmon, S. B. Libicki, and C. R. Robertson, “A theoretical investigation of convective transport in the hollow-fiber reactor,” Chemical Engineering Communications, vol. 66, no. 1, pp. 221–248, 1988.

J. Lu and W. Q. Lu, “A numerical simulation for mass transfer through the porous membrane of parallel straight channels,” International Journal of Heat and Mass Transfer, vol. 53, no. 11-12, pp. 2404–2413, 2010.

J. L. Anderson, “Configurational effect on the reflection coefficient for rigid solutes in capillary pores,” Journal of Theoretical Biology, vol. 90, no. 3, pp. 405–426, 1981.

K. Annan, “Mathematical modeling of the dynamic exchange of solutes during bicarbonate dialysis,” Mathematical and Computer Modelling, vol. 55, no. 5-6, pp. 1691–1704, 2012.

Y. Lee, R. D. Noble, B. Y. Yeom, Y. I. Park, and K. H. Lee, “Analysis of removal by hollow fiber membrane contactors,” Journal of Membrane Science, vol. 194, no. 1, pp. 57–67, 2001.

P. V. Dankwerts and M. M. Sharma, “The absorption of carbon dioxide into solutions of alkalis and amines (with some notes on hydrogen sulfide and carbonyl sulfide),” Chemical engineering, 244, 1966

K. Soudan, E. S. Ricanati, J. B. Leon, and A. R. Sehgal, “Determinants of metabolic acidosis among hemodialysis patients,” Hemodialysis International, vol. 10, no. 2, pp. 209–214, 2006.

S. R. Suchdeo and J. S. Schultz, “The permeability of gases through reacting solutions: the carbon dioxide-bicarbonate membrane system,” Chemical Engineering Science, vol. 29, no. 1, pp. 13–23, 1974.

J. F. Patzer, “Oxygen consumption in a hollow fiber bioartificial liver-revisited,” Artificial Organs, vol. 28, no. 1, pp. 83–98, 2004.

M. Labecki, J. M. Piret, and B. D. Bowen, “Two-dimensional analysis of fluid flow in hollow-fibre modules,” Chemical Engineering Science, vol. 50, no. 21, pp. 3369–3384, 1995.

S. Tong and F. Yuan, “An equivalent length model of microdialysis sampling,” Journal of Pharmaceutical and Biomedical Analysis, vol. 28, no. 2, pp. 269–278, 2002.

M. Külz, B. Nederlof, and H. Schneider, “In vitro and in vivo evaluation of a new dialyzer,” Nephrology Dialysis Transplantation, vol. 17, no. 8, pp. 1475–1479, 2002.

W. Ding, L. He, G. Zhao, H. Zhang, Z. Shu, and D. Gao, “Double porous media model for mass transfer of hemodialyzers,” International Journal of Heat and Mass Transfer, vol. 47, no. 22, pp. 4849–4855, 2004.

J. L. Anderson, “Configurational effect on the reflection coefficient for rigid solutes in capillary pores,” Journal of Theoretical Biology, vol. 90, no. 3, pp. 405–426, 1981.

F. J. Gennari, “Acid-base considerations in end stage renal disease,” in Principles and Practice of Dialysis, W. L. Henrich, Ed., pp. 393–407, Lippincott Williams & Wilkins, Philadephia, Pa, USA, 2004.