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

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, HCO₃⁻ concentration gradients peaked at the same position. Also, CO₂ concentration decreased rapidly within the first 47 minutes while optimal HCO₃⁻ 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 𝜕𝜌𝑐^UnsteadyTerm+𝜕𝑡^{ConvectionTerm}div𝜌𝑐𝐮=div𝐷⋅grad^{DiffusionTerm}(𝑐)+𝑆^SourceTerm_𝑐, (1) 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 div() and grad() 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 (div{𝜌𝑐𝐮}) plus the rate they are being added by diffusion (div{𝐷grad(𝑐)}) 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 𝑆_𝑐=0.(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 grad(𝑐) except for the radial direction. Using assumptions 1–4, (1) becomes 𝜌⋅div(𝑐𝐮)=𝐷⋅divgrad(𝑐). (2)Integrating (2) and transforming the resulting equation using the divergence theorem gives 𝜌_𝐶𝑉div(𝑐𝐮)𝑑𝑉=𝐷_𝐶𝑉(𝑐)𝜌divgrad𝑑𝑉,_𝐹(𝐧⋅𝑐𝐮)𝑑𝐴=𝐷_𝐹𝐧⋅grad(𝑐)𝑑𝐴, (3) 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 𝜌𝑢_𝑟𝑐𝑟+𝛿𝑟2−𝑐𝑟−𝛿𝑟2=𝐷𝜕𝑐𝜕𝑟𝑟+𝛿𝑟2−𝜕𝑐𝜕𝑟𝑟−𝛿𝑟2. (4) 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) 𝜌𝑢_𝑟𝑐𝑟+𝛿𝑟2−𝐷𝜕𝑐𝜕𝑟𝑟+𝛿𝑟2=𝜌𝑢_𝑟𝑐𝑟−𝛿𝑟2−𝐷𝜕𝑐𝜕𝑟𝑟−𝛿𝑟2. (5) 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 𝜌𝑢_𝑟𝑐(𝑟)−𝐷𝜕𝑐(𝑟)𝜕𝑟=𝐽_𝑠. (6) The exact solution for (6) in terms of the concentration 𝑐(𝑟) is of the form 𝑐(𝑟)=𝑘₁𝑒^𝑘2𝑟+𝑘₃. (7) Substituting (7) into (6) and gathering terms 𝜌𝑢_𝑟𝑐(𝑟)−𝐷𝜕𝑐(𝑟)𝜕𝑟=𝜌𝑢_𝑟𝑘₁𝑒^𝑘2𝑟+𝜌𝑢_𝑟𝑘₃−𝐷𝑘₁𝑘₂𝑒^𝑘2𝑟=𝑘₁𝜌𝑢_𝑟−𝐷𝑘₂𝑒^𝑘2𝑟+𝜌𝑢_𝑟𝑘₃=𝐽_𝑠. (8) The coefficient of the exponential must be zero and the constant term equal 𝐽_𝑠, thus 𝑘₂=𝜌𝑢_𝑟𝐷,𝑘₃=𝐽_𝑠𝜌𝑢_𝑟. (9) 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 𝑐(𝑟)=𝛼₃−𝛼₁𝑒^−𝛼2𝑟, (10) 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 (10⁻⁴) fiber pore size [30], we describe the flow passing through the membrane (see Figure 1(b)) by simplified Kedem-Katchalsky (K-K) equations 𝐽_𝑣≈𝐿_𝑝𝐽Δ𝑃,_𝑠≈𝐶^∗_𝑠𝐽_𝑣+𝑃_𝑠Δ𝑐_𝑠. (11)𝐽_𝑣 (m/s) is ultrafiltration velocity or volumetric flux across the membrane; 𝐽_𝑠(kg/m²s) is solute flux across the membrane; 𝐿_𝑝(m/sPa) is the hydraulic permeability of the membrane; 𝑃_𝑠(m/s) is solute diffusive permeability coefficient of a membrane; 𝑐^∗_𝑠(kg/m³) represents the average solute concentration at each side of the membrane; Δ𝑐_𝑠(kg/m³) is solute concentration difference (i.e., transmembrane concentration) across the membrane. The parameter Δ𝑃(Pa) is the membrane surface hydraulic permeability of the membrane. Thus, the membrane interfacial conditions for the blood-side model are 𝑢_𝑧=0,𝑢_𝑟=𝐽_𝑣,𝐷_𝑠𝜕𝑐_𝑠𝜕𝑟=𝐽_𝑣𝑐_𝑠−𝐽_𝑠. (12)

(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., 0≤𝑧≤𝐿) and 𝑟 is taken along the radial direction. An axisymmetric domain, where 𝑟 is chosen to lie in the range 0<𝑟<𝑟_𝑏 between 𝑧=0 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] 1𝑟𝜕𝑟𝑢_𝑟+𝜕𝑟𝜕𝑢_𝑧𝑢𝜕𝑧=0,_𝑟𝜕𝑢_𝑟𝜕𝑟+𝑢_𝑧𝜕𝑢_𝑟1𝜕𝑧=−𝜌𝜕𝑝+𝜇𝜕𝑟𝜌1𝑟𝜕𝑟𝜕𝑟𝜕𝑢_𝑟−𝑢𝜕𝑟𝑟^𝑟2+𝜕²𝑢_𝑟𝜕𝑧²,𝑢_𝑟𝜕𝑢_𝑧𝜕𝑟+𝑢_𝑧𝜕𝑢_𝑧1𝜕𝑧=−𝜌𝜕𝑝+𝜇𝜕𝑧𝜌1𝑟𝜕𝑟𝜕𝑟𝜕𝑢_𝑧+𝜕𝜕𝑟²𝑢_𝑧𝜕𝑧², (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 𝑧=0 and 0<𝑟<𝑟_𝑏 are obtained [31, 34] 𝑢_𝑟(𝑟)=0,𝑢_𝑧(𝑟)=2𝑞_𝑏𝑁𝜋𝑟⁴_𝑏𝑟²_𝑏−𝑟². (14) 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 𝑟=0𝑢_𝑟=𝑢_𝑧=0,𝑢_𝑟=𝜕𝑢_𝑧𝜕𝑟=0at𝑟=0;0≤𝑧≤𝐿. (15)

The convection-diffusion equation governing the mass transport of solutes 𝑠 coupled to the blood velocity field is given by 𝜕𝑐^{TransientTerm}_𝑠+𝜕𝑡^{ConvectiveTerm}𝑢_𝑧𝜕𝑐_𝑠𝜕𝑧+𝑢_𝑟𝜕𝑐_𝑠=𝜕𝑟_{DiffusiveTerm}𝐷_𝑠𝜕²𝑐_𝑠𝜕𝑟²+1𝑟𝜕𝑐_𝑠+𝜕𝜕𝑟²𝑐_𝑠𝜕𝑧²+𝐵^BufferTerm_𝑠, (16) 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 𝑐_𝑠(𝑧,𝑟,0)=𝑐_𝑠0,𝑐_𝑠(0,𝑟,𝑡)=𝑐_𝑠0,𝜕𝑐_𝑠(𝑧,0,𝑡)𝜕𝑟=0. (17)𝐵_𝑠 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 𝑠=[CO₂,HCO₃⁻] given by [31, 35] 𝐵_CO2=−𝑘₊21+𝛼CO₃²⁻HCO₃⁻CO₂−𝛽HCO₃⁻CO₃²⁻,𝐵_HCO3⁻=−2𝐵_CO2, (18) where 𝛼=𝑘₂𝐾₄/2𝑘₊𝐾₃,  𝛽=𝑘₋𝐾₃/𝑘₊,   and  [HCO₃⁻]/[CO₃⁻²] = 20. The rate controlling reactions for the carbonate and bicarbonate ions are given as [31, 35] CO₂+H₂O⇌^𝑘+_𝑘−H⁺+HCO₃⁻,CO₂+OH⇌^−𝑘2_𝑘−2HCO₃⁻, (19) where 𝑘₊,𝑘₋,𝑘₂, and 𝑘₋₂ are their reaction constants with their equilibrium constants defined as 𝐾₁ and 𝐾₂, respectively. The overall reaction is CO₂+CO₃²⁻+H₂O⇌2HCO₃⁻, (20) with the following fast reactions assumed to be at equilibrium HCO⇌₃^−𝐾3H⁺+CO₃⁻,H₂O+⇌^𝐾4H⁺+OH⁻, (21) 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.

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 𝑢_𝑧=0 at 𝑟=0 and 𝑟=𝑟_𝑏𝑢_𝑟𝑟_𝑏𝑞=−^𝑑2𝜋𝑟_𝑏𝐿_𝑟,∀𝑟_𝑏,𝑢_𝑧=2𝑞_𝑑𝜋𝑟²_𝑑𝐿_𝑟ln𝑟/𝑟_𝑑−𝑟/𝑟_𝑑²−1/𝜅²−1ln(𝜅)𝜅²+1ln(𝜅)+1−𝜅². (22) 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] 𝑅_𝑠=𝜀𝑐_𝑠𝑐_𝑠0−𝑐_𝑠, (23) 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 𝜕𝑐^Transient_𝑠+𝜕𝑡^Convective𝑢_𝑧𝜕𝑐_𝑠𝜕𝑧+𝑢_𝑟𝜕𝑐_𝑠=𝜕𝑟_Diffusive𝐷_𝑠𝜕²𝑐_𝑠𝜕𝑟²+1𝑟𝜕𝑐_𝑠+𝜕𝜕𝑟²𝑐_𝑠𝜕𝑧²+𝑅^{Replenishment}_𝑠. (24)

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 𝜕𝜕𝑡_𝐶𝑉𝑐_𝑠𝑑𝑉+_𝐶𝑉𝑢_𝑧𝜕𝑐_𝑠=𝜕𝑧𝑑𝑉_𝐶𝑉𝐷_𝑠𝜕²𝑐_𝑠𝜕𝑟²+1𝑟𝜕𝑐_𝑠+𝜕𝜕𝑟²𝑐_𝑠𝜕𝑧²𝑑𝑉+_𝐶𝑉𝜓𝑑𝑉. (25) 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 𝑐_𝑠0𝑐𝑢=_𝑠𝑢_𝑧(𝑟). (26) Therefore, the divergence of the convective flow vector 𝑐_𝑠𝑢 using cylindrical coordinate is div𝑐_𝑠𝑢=1𝑟𝜕𝑟𝑐_𝑠𝑢_𝑟+1𝜕𝑟𝑟𝜕𝑐_𝑠𝑢_𝜃+𝜕𝜃𝜕𝑐_𝑠𝑢_𝑧=𝜕𝑧𝜕𝑐_𝑠𝑢_𝑧=𝜕𝑧𝜕𝑐_𝑠𝑢𝜕𝑧^𝑧+𝑐_𝑠𝜕𝑢_𝑧=𝜕𝑧𝜕𝑐_𝑠𝑢𝜕𝑧^𝑧, (27) 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 𝜕𝑢_𝑧/𝜕𝑧=0. Thus, the second term on the left-hand side of (25) is _𝐶𝑉𝑢_𝑧𝜕𝑐_𝑠𝜕𝑧𝑑𝑉=_𝐶𝑉div𝑐_𝑠𝑢𝑑𝑉=_𝐴𝑛⋅𝑐_𝑠𝑢𝑑𝐴. (28) 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 𝑛⋅𝑐_𝐴𝑠𝑢𝑑𝐴=_+𝑧𝑐_𝑠𝑟,𝑧_𝑐+Δ𝑧2𝑢_𝑧(𝑟)−𝑟𝑑𝑟𝑑𝜃_−𝑧𝑐_𝑠𝑟,𝑧_𝑐−Δ𝑧2𝑢_𝑧(𝑟)𝑟𝑑𝑟𝑑𝜃, (29) 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, 𝑛⋅𝑐_𝐴𝑠𝑢𝑑𝐴=𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧2𝑢_𝑧𝑎(𝑟)_+𝑧𝑟𝑑𝑟𝑑𝜃−𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧2𝑢_𝑧𝑎(𝑟)_−𝑧𝑟𝑑𝑟𝑑𝜃,=𝑢_𝑧𝑎(𝑟)𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧2−𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧2Δ𝐴_𝑧. (30)

3.1.2. Diffusive Term

Since diffusion is driven by concentration gradient, the divergence vector field in cylindrical coordinates is (𝑐)=1divgrad𝑟𝜕𝑟𝜕𝑟𝜕𝑐+𝜕𝜕𝑟𝜕𝑧𝜕𝑐,=1𝜕𝑧𝑟1⋅𝜕𝑐𝜕𝜕𝑟+𝑟⋅²𝑐𝜕𝑟²+𝜕²𝑐𝜕𝑧²,=1𝑟𝜕𝑐+𝜕𝜕𝑟²𝑐𝜕𝑟²+𝜕²𝑐𝜕𝑧². (31) Using the divergence theorem, the diffusion term in (25) could be written as 𝐷_𝑠𝐶𝑉divgrad𝑐_𝑠𝑑𝑉=𝐷_𝑠𝐴𝑛⋅𝜕𝑐_𝑠𝜕𝑟𝜕𝑐_𝑠𝜕𝑧𝑑𝐴. (32) 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 𝐷_𝑠𝐴𝑛⋅𝜕𝑐_𝑠𝜕𝑟𝜕𝑐_𝑠𝜕𝑧𝑑𝐴=𝐷_𝑠+𝑟𝜕𝑐_𝑠𝑟_𝑐+Δ𝑟/2,𝑧𝜕𝑟𝑑𝐴_𝑟−_−𝑟𝜕𝑐_𝑠𝑟_𝑐−Δ𝑟/2,𝑧𝜕𝑟𝑑𝐴_𝑟+_+𝑧𝜕𝑐_𝑠𝑟,𝑧_𝑐+Δ𝑧/2𝜕𝑧𝑑𝐴_𝑧−_−𝑧𝜕𝑐_𝑠𝑟,𝑧_𝑐−Δ𝑧/2𝜕𝑧𝑑𝐴_𝑧. (33) 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 𝐷_𝑠𝐴𝑛⋅𝜕𝑐_𝑠𝜕𝑟𝜕𝑐_𝑠𝜕𝑧𝑑𝐴=𝐷_𝑠Δ𝐴_𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐+Δ𝑟/2,𝑧−𝜕𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐−Δ𝑟/2,𝑧𝜕𝑟+Δ𝐴_𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧/2−𝜕𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧/2.𝜕𝑧 (34) Thus, the LHS of (25) with 𝑐_𝑠 defining the volume average of solute concentration is 1Δ𝑉_𝐶𝑉𝜕𝑐_𝑠+1𝜕𝑡𝑑𝑉Δ𝑉_𝐶𝑉𝑢_𝑧𝜕𝑐_𝑠𝜕𝜕𝑧𝑑𝑉=1𝜕𝑡Δ𝑉_𝐶𝑉𝑐_𝑠+𝑢𝑑𝑉_𝑧𝑎(𝑟)Δ𝐴_𝑧𝑐Δ𝑉^𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧2−𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧2=𝜕𝑐_𝑠𝜕𝑡+𝑢_𝑧𝑎(𝑟)𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧/2−𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧/2.Δ𝑧 (35) The RHS of (25) becomes 1Δ𝑉_𝐶𝑉𝐷_𝑠𝜕²𝑐_𝑠𝜕𝑟²+1𝑟𝜕𝑐_𝑠+𝜕𝜕𝑟²𝑐_𝑠𝜕𝑧²1𝑑𝑉+Δ𝑉_𝐶𝑉=𝐷𝜓𝑑𝑉_𝑠×Δ𝑉Δ𝐴_𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐+Δ𝑟/2,𝑧−𝜕𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐−Δ𝑟/2,𝑧𝜕𝑟+Δ𝐴_𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧/2−𝜕𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧/2+𝜕𝑧𝜓,=𝐷_𝑠1Δ𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐+Δ𝑟/2,𝑧−𝜕𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐−Δ𝑟/2,𝑧+1𝜕𝑟Δ𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧/2−𝜕𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧/2+𝜕𝑧𝜓, (36) where 𝜓=(1/Δ𝑉)∫_𝐶𝑉𝜓𝑑𝑉.

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𝑟^∗=𝑟𝐿;𝑧^∗=𝑧𝐿;𝑢^∗_𝑟=𝑢_𝑟𝑈;𝑢^∗_𝑧=𝑢_𝑧𝑈;𝐴₁=𝑘₊𝐿𝑈;Sh_𝑠=𝐿𝑃_𝑠𝐷_𝑠;𝑐^∗_𝑠=𝑐_𝑠𝑐_𝑠0;𝑡^∗=𝑡𝑈𝐿𝑟;𝜙=^𝑏𝐿;Pe=𝐿𝑈𝐷_𝑠;Re=𝜌𝑈𝑟_𝑏𝜇;𝐸_𝑠=𝐿_𝑝𝐿Δ𝑝2𝐷_𝑠; (37) where Pe and Sh are Pe’clet and Sherwood numbers, respectively, and the ratio of momentum diffusivity and mass diffusivity is donated by 𝐸. Pe=𝐿𝑈/𝐷_𝑠=𝑞_𝑏/𝜋𝜙𝑟_𝑏𝐷_𝑠 expressed the relative importance of convection to diffusion while Re=𝜌𝑈𝑟_𝑏/𝜇=𝜌𝑞_𝑏/𝜋𝜇𝑟_𝑏 related inertial effects to viscous effects. Since dialysis devices employ laminar fluids flow with Re≪1 the inertial effects would be irrelevant [31, 41].

Table 2 

The reference variables with their description.

Substituting the dimensionless variables in (37) into (35) and (36), simplifying notations, and dropping the superscript “*” resulted in 𝜕𝑐_𝑠𝜕𝑡+𝑢_𝑧𝑎(𝑟)𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧/2−𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧/2=1Δ𝑧1PeΔ𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐+Δ𝑟/2,𝑧−𝜕𝑟𝜕𝑐_𝑠𝑎𝑟_𝑐−Δ𝑟/2,𝑧+1𝜕𝑟Δ𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧/2−𝜕𝑧𝜕𝑐_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧/2+𝐿𝜕𝑧𝑈𝑐_𝑠0𝜓(𝑟,𝑧). (38) The dimensionless initial and boundary conditions, buffer and replenishment, and membrane interfacial conditions corresponding to (38) were as follows.

Blood-Side Inlet Velocity Conditions𝑢_𝑟=0,𝑢_𝑧=2𝑁𝑟1−²𝜙². (39)

Dialysate-Side Inlet and Outlet Velocities𝑢_𝑟𝑟=−^𝑑2𝜅𝐿_𝑟∀𝑟_𝑏𝑢_𝑧=2𝐿_𝑟𝜙²𝜅²−1ln(𝜅)⋅ln(𝜅𝑟)−ln𝜙−(𝜅𝑟)²−𝜙²𝜙²𝜅⁴−1ln(𝜅)−𝜅²−1²ln(𝜅) (40)

No Slip and Axisymmetric Conditions𝑢_𝑟=𝑢_𝑧=0,𝑢_𝑟=𝜕𝑢_𝑧𝜕𝑟=0at𝑟=0;0≤𝑧≤1. (41)

Inlet and Outlet Blood and Dialysate Concentrations𝑐_𝑠(𝑧,𝑟,0)=1,𝑐_𝑠(0,𝑟,𝑡)=1,𝜕𝑐_𝑠(𝑧,0,𝑡)𝜕𝑟=0. (42)

Buffer and Replenishment Terms𝜓_CO2𝐴=−¹𝑐_𝑠0CO₂−20𝛽(1+0.1𝛼),𝜓_HCO3−=−2𝜓_CO2,𝜓_𝑅𝑠=𝜀𝑐_𝑠0𝐿𝑈𝑐_𝑠1−𝑐_𝑠, (43) where 𝜓_CO2and𝜓_HCO3− represented dimensionless buffer terms in blood side and 𝜓_𝑅𝑠 depicted dimensionless replenishment term for solute 𝑠=CO₂ and HCO₃⁻.

Blood-Membrane Interfacial Conditions𝑢_𝑧=0,𝑢_𝑟=𝐽_𝑣𝑈=𝐸_𝑠Pe_𝑠,𝜕𝑐_𝑠𝜕𝑟=𝐸_𝑠1−𝐶^∗_𝑠−Sh_𝑠⋅Δ𝑐_𝑠, (44) where Sh_𝑠 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., 6.65  ×  10⁻⁶ m³ s⁻¹) and 800 mL/min (1.33×10⁻⁵ 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.

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 𝑅_𝑏max, 𝑅_𝑑max, and 𝑅_max. 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 𝑧=0, so the first internal grid point is at 0+𝑧_grid/2, where 𝑧_grid 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 𝑧=−𝑧_grid/2. Therefore, the grid size in each axis was 𝑧_grid=1𝑧_max−2,𝑅_𝑏grid=𝑅_𝑏𝑅_𝑎𝑅_𝑏max,𝑅−2_𝑑grid=𝑅_𝑑𝑅_𝑎𝑅_𝑑max,−2 (45) since the domain of interest in the model was 0<𝑧<1 for the 𝑧-direction and 0<𝑟<𝑟_𝑏/𝐿 and 0<𝑟<𝑟_𝑑/𝐿 in the 𝑟-direction for blood and dialysate sides, respectively.

The indices of the variables (such as solute concentration) ran from 𝑖=1 to 𝑧_max and 𝑗=1 to 𝑅_𝑏max or 𝑅_𝑑max. Therefore, the boundaries of the spatial domains were between indices 𝑗=1 and 2 and between 𝑗=𝑅_𝑏max−1 and 𝑅_𝑏max (blood side), similarly for the dialysate side, it was between 𝑖=1 and 2 and between 𝑖=𝑧_max−1 and 𝑖=𝑧_max. Since the boundary conditions define the values of the variables there, the numerical model only calculated the values in the range from 2 to 𝑧_max−1 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 𝑡(𝑘)=(𝑘−1)𝑑𝑡,𝑘=1,2,…. 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.

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: 𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧2=𝑐_𝑏𝑟,𝑧_𝑐=𝑐_𝑏,𝑐𝑖,𝑗,𝑘+1,𝑠_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧2=𝑐_𝑏𝑟,𝑧_𝑐−Δ𝑧=𝑐_𝑏𝑖−1,𝑗,𝑘+1,𝑠. (46)

Dialysate Side: 𝑐_𝑠𝑎𝑟,𝑧_𝑐+Δ𝑧2=𝑐_𝑑𝑟,𝑧_𝑐+Δ𝑧=𝑐_𝑑,𝑐𝑖+1,𝑗,𝑘+1,𝑠_𝑠𝑎𝑟,𝑧_𝑐−Δ𝑧2=𝑐_𝑑𝑟,𝑧_𝑐=𝑐_𝑑𝑖,𝑗,𝑘+1,𝑠. (47)

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 𝑡=70 minutes.

(a) CO2 Concentration gradients at various times

(b) HCO3− 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 𝑡=60 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 (𝑧=20 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 CO₂ and HCO₃⁻ 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 CO₂ may be caused by (i) the active hydrogen ion from the blood reacting with bicarbonate ions to produce more of the CO₂ and/or (ii) the incomplete dissociation of bicarbonate ion into CO₂ 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) CO2 Diffusion fluxes at various radial positions

(b) HCO3− Diffusive fluxes at various radial positions

Figure 5 

Bicarbonate and carbon dioxide diffusive fluxes in the membrane at dialyzer axial distance 𝑧=20 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 (CO₂ and HCO₃⁻) during a high-flux hollow-fiber dialysis session. The magnitude of CO₂ and HCO₃⁻ 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) CO2 convective fluxes at various radial positions

(b) HCO3− 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) CO2 Convective and total fluxes in membrane

(b) HCO3− 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 CO₂ concentration profiles at different time periods were shown as a function of membrane axial distance. The result clearly showed that the CO₂ concentration decreased as the axial length increased and that most CO₂ desorption occurred within 47 minutes of the therapy session. The constant profiles after 47 minutes indicated that the membrane was fully saturated with CO₂ concentration during the therapy. In Figure 8(b), CO₂ concentration profiles as a function of the first 70 minutes of dialysis session was shown at different radial positions. The result demonstrated that CO₂ concentration increased as one moved toward the end of the membrane. Also, the CO₂ 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 CO₂ saturation in the membrane.

(a) CO2 Concentrations at different time periods

(b) CO2 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 HCO₃⁻ concentration increased as the membrane axial distance increased (see Figure 9(a)). Also, it was shown that HCO₃⁻ concentration increased rapidly within the first 30 minutes and then became stable after 40 minutes of the therapy session. In Figure 9(b), the HCO₃⁻ concentration profiles at axial distance (𝑧=0.20m) against time for various radial distances was shown. HCO₃⁻ concentration increased and was nearly saturated in the membrane. The rate of increased and the degree at which HCO₃⁻ 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 HCO₃⁻ concentration levels during dialysis therapy.

(a) HCO3− Concentrations at various time periods

(b) HCO3− 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 CO₂ concentration gradients appeared to diminish from their maxima and shift toward the end of the membrane, HCO₃⁻ concentration gradients peaked at the same position. In addition, the magnitude of the HCO₃⁻ concentration gradient was large enough to activate HCO₃⁻ 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 (CO₂ and HCO₃⁻) 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 CO₂ within 47 minutes. Further investigation showed an increased CO₂ concentration towards the end of the membrane which indicated higher CO₂ concentration near the walls of the fiber than the fiber center at the same axial distance. Similarly, there was an optimal HCO₃⁻ 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.