Investigating the Dialysis Treatment Using Hollow Fiber Membrane: A New Approach by CFD

Due to the increase in the number of people affected by chronic renal failure, the demand for hemodialysis treatment has increased considerably over the years. In this sense, theoretical and experimental studies to improve the equipment (hemodialyzer) are extremely important, due to their potential impact on the patient’s life quality undergoing treatment. To contribute to this research line, this work aims to study the fluid behavior inside a hollow fiber dialyzer using computational fluid dynamics. In that new approach, the blood is considered as multiphase fluid and the membrane as an extra flow resistance in the porous region (momentum sink). The numerical study of the hemodialysis process was based on the development of a mathematical model that allowed analyzing the performance of the system using Ansys® Fluent software. The predicted results were compared with results reported in the literature and a good concordance was obtained. The simulation results showed that the proposed model can predict the fluid behavior inside the hollow fiber membrane adequately. In addition, it was found that the clearance decreases with increasing radial viscous resistance, with greater permeations in the vicinity of the lumen inlet region, as well as the emergence of the retrofiltration phenomenon, characteristic of this type of process. Herein, velocity, pressure, and volumetric fraction fields are presented and analyzed.


Introduction
In Brazil, according to the Brazilian nephrology census carried out in July 2016, approximately 122,825 people are on dialysis treatment due to renal insufficiency. Unevenly distributed throughout the country, there are 747 active dialysis units (4% in the north region, 18% in the northeast region, 7% in the midwest, 49% in the southeast, and 22% located in the south region), responsible for providing this treatment for the population [1].
Renal insufficiency can be defined as the loss of the kidney's ability to maintain the electrolyte balance of the body and remove the metabolic breakdown products. In its acute phase, the disease causes a rapid reduction in renal function, being reversible if treated properly. However, chronic renal insufficiency, unlike acute, is characterized by the gradual and irreversible loss of these functions [2].
The progressive accumulation of compounds through metabolism, normally excreted by healthy kidneys, is usually known as the uremic syndrome. This symptom of renal insufficiency can lead to a reduction in the filtration rate of the kidneys to values below 5 mL/min [3,4]. The incidence of this disease has gradually increased in recent years, as a reflection of the population aging that needs to live each day with several comorbidities [5].
Patients suffering from renal insufficiency can be treated initially through blood pressure control, medications, and dietary treatments. However, over time, dialysis treatment or kidney transplantation is necessary [2,6]. In hemodialysis therapy, a device known as a dialyzer is used to supply the kidneys' function. This equipment consists of a hull (shell) fitted with a bundle of tubes (hollow fiber membranes). The dialyzer allows blood to flow through the tubes, enabling the removal of metabolic waste by mass transfer (diffusion and convection) that occurs through a porous membrane between the blood and the dialysis solution (dialysate), contained on the shell side [7].
Clark et al. [8] mention that the dialyzer is equipment consisting basically of three components: the dialysate fluid compartment, the blood compartment, and the membrane. The authors emphasize that the knowledge about the physicochemical phenomena and mechanisms involved in this process contribute to optimizing the operational parameters involved in hemodialysis treatment.
Ding et al. [9], when studying an artificial kidney, observed that the simulation of a dialyzer based on the finite-element method is a viable alternative to investigate the behavior of the concentration and velocity fields inside the hollow fiber membrane, enabling the obtainment of expressive results, useful for predicting the removal rate of toxins from the blood.
Some works have been reported in the literature evaluating the phenomenon of mass transfer in dialyzers (hollow fiber membranes). Among the various experimental works, we can mention Klein et al. [10] and Liao et al. [11]. Among the theoretical works, we can mention Lu and Lu [12], which proposed to analyze the mass transfer in a dialysis system with the concurrent flow in parallel plates. Among the numerical works, we can mention Pstras et al. [13], who presented the main mathematical models used to study/optimize hemodialysis therapy, and Cancillaet al. [14], who developed a CFD model to simulate the hemodialysis process in hollow fiber membrane modules.
In addition to the aforementioned research, numerical works can be found discussing dialysis treatment from the perspective of chemical species transfer, such as Gostoli and Gatta [15], who studied mass transfer in countercurrent and concurrent flow in a capillary; Ding et al. [16], when developing a double porous zone model for mass transfer in a hemodialyzer, and Kanchan and Maniyeri [17], Liao et al. [18], Lu and Junfeng Lu, [6] and Donato et al. [19] when studying the solute transport in hemodialysis membranes with a two-dimensional approach. However, to investigate the dialysis process from the perspective of momentum transport, Choi et al. [20] and Kim et al. [21] modeled the blood as a multiphase fluid.
Hemodialysis is the main treatment for patients with chronic renal insufficiency. The procedure is performed in specialized nephrology services and has an average duration of three to four hours, requiring three weekly dialysis sessions, which makes the treatment very difficult for human beings. Therefore, contributions in this area represent an expectation to improve the quality of life of patients with this comorbidity. Although works are being developed in this area, investigations about fluid dynamics inside the hollow fiber are still scarce, especially dealing with blood flow through a multiphase approach. Therefore, in addition to the aforementioned research, this work aims to investigate the hemodialysis process in a hollow fiber membrane via computational fluid dynamics.
Hemodialysis is the main treatment for patients with chronic renal insufficiency. The procedure is performed in specialized nephrology services and has an average duration of three to four hours, requiring three weekly dialysis sessions, which makes the treatment very difficult for human beings. Therefore, contributions in this area represent an expectation to improve the quality of life of patients with this comorbidity. Although works are being developed in this area, investigations about fluid dynamics inside the hollow fiber are still scarce, especially dealing with blood flow through a multiphase approach, where the blood is modeled as two distinct phases (blood and contaminant). Therefore, in addition to the aforementioned research, this work aims to investigate the hemodialysis process in a hollow fiber membrane via computational fluid dynamics. The innovative aspect in this research deals with each phase's mathematical consideration. Herein, each involved phase was defined as an identifiable class of material that has a particular inertial response and interaction with the flow field in which it is immersed. Besides, the porous medium is considered as an extra resistance to flow in the porous region in the form of a moment sink by adding a source term to the momentum equation applied in each element in that region. Clearly, it is an attempt to predict the hemodialysis process based in a new mathematical procedures and tools. Thus, the authors strongly recommend that more sophisticated approaches must be investigated from this research.

Problem Description
The research is based on a model CT190G dialyzer (Baxter Healthcare Co., McGaw Park, IL, USA) consisting of a hull and a bundle of 12,000 triacetate cellulose hollow fiber. A schematic of the equipment and details of the membrane section used in the methodology can be seen in Figure 1, as described by Liao et al. [18]. The hollow fiber membrane consists of three parts: the shell (dialysate flow region), the porous membrane, and the lumen (blood flow region). Inside the equipment, blood (water and urea contaminant) flows into the lumen domain and the dialysate into the shell domain, with the flow in opposite directions (countercurrent), as illustrated in Figure 2. The dimensions of the equipment are listed in Table 1. Table 1. Dimensions of the hollow fiber membrane [18].

Computational Domain
To perform the numerical simulations, three numerical meshes (M 1 , M 2 , and M 3 ) were built using the Ansys ® Designer Modeler and Meshing 15.0 software (Canonsburg, PA, USA), as illustrated in Figure 3.

Mathematical Modeling
For the study of the hemodialysis process using a hollow fiber membrane section, the following considerations were assumed: The proteins present in the blood were disregarded; • Adsorption of urea on the membrane contact surface, blockage of membrane pores, formation of concentration polarization layer, and chemical reactions are disregarded; • Only one section of the hollow fiber membrane is considered, due to the angular symmetry presented by the geometry; • The Eulerian-Eulerian approach was adopted for multiphase flow.
After due consideration, the mass conservation equation for phase q (Equation (1)) and linear momentum (Equation (2)) can be written as:

•
Mass conservation equation for the non-porous media where f is the volume fraction of the phase q, → ν q is the velocity vector of the phase q, and ρ is the phase density.

•
Linear momentum equation where P is the pressure shared by all phases and τ q is the stress-strain tensor. The stress-strain tensor is defined by Equation (3) as follows: where µ q and λ q are the viscosity and shear of phase q.
The term referring to the interface forces, → R pq , depends on pressure, friction, cohesion, and other acting effects, and is subject to the conditions of Equations (4) and (5).
Ansys Fluent ® software uses an interaction term between the forces, described by: where K pq is the interface moment exchange term, given by: where f a is the drag function, and τ p is the particulate relaxation time, defined as: where C D is the drag coefficient, and Re is the relative Reynolds number, defined for the primary phase (q) and for the second phase (p), and d p is the diameter of the droplet or bubble of the secondary phase. The drag coefficient and the relative Reynolds number are calculated using Equations (11) and (12) [22].
• Linear momentum equation for the porous medium The porous media model is formed by incorporating extra flow resistance in the porous region, in the form of a momentum sink. This occurs by adding a source term in the momentum equation, applied to the elements of that region. The added source term is composed of two parts: a viscous loss term (Darcy) and an inertial loss term, as seen in Equation (6).
where S i is the source term for the i (x, y, or z) momentum equations, |ν| is the magnitude of the velocity, and µ is the viscosity. The source term added (momentum sink) contributes to the pressure gradient in the porous cell, creating a pressure drop that is proportional to the fluid velocity. For simple homogeneous porous media, Equation (13) can be rewritten as: where α is the permeability and C 2 is the inertial resistance factor. When using the porous medium model, it is considered that the porous cells are completely open, and the only resistance imposed to the flow is given in the form of the viscous (1/α) and inertial (C 2 ) resistance coefficients.
For laminar flow through porous media, the constant C 2 can be considered zero, as the pressure drop is normally proportional to velocity. By disregarding diffusion and convective acceleration, the porous medium model can be expressed by Darcy's Law, as follows: Therefore, the pressure gradient inside the porous region can be given by: where ∆n x , ∆n y , and ∆n z are the thicknesses of the porous medium in the x, y, and z directions, respectively.

Conditions Used in Simulations (a) Initial and boundary conditions
To complete the mathematical modeling, the following initial and boundary conditions illustrated in Figure 4 were used: • Initial conditions The initial concentration of the contaminant (urea) at the entrance of the lumen, C in , at t = 0 s, so that: where C in = 0.7 kg/m 3 . In addition, the fiber was considered to be completely filled with saline solution (water) at time t 0 , before starting the process.

• Boundary conditions
On the axis of symmetry The following conditions were admitted on the axis of symmetry: In the domain inlets A condition of constant values for blood and dialysate volumetric fluxes at the domain inlets, in countercurrent, was assumed.
where Q Bin is the blood volumetric flowrate (water + urea), Q Din is the dialysate volumetric flowrate, and N is the number of dialyzer tubes proposed by Liao et al. [18].
In the domain outlets A zero-pressure condition at the mass flow outlets was assumed.
(b) Thermophysical parameters of membrane and fluids The thermo-physical parameters of the fluids and membrane are shown in Table 2. Table 2. Thermo-physical properties and parameters of fluids and membrane [18].

Studied Cases
To evaluate the influence of the mesh on the simulation results, the grid convergence index (GCI) method was applied to three meshes with different element densities, as shown in Table 3. In the simulations performed, the lumen feed flux, Q Bin , the shell feed flux, Q Din , the axial viscous resistance, 1/α x , the radial viscous resistance, 1/α y , and the urea concentration in the lumen feed, C in , were kept constant (Table 4). After the mesh study was completed, the optimized mesh was selected. With the appropriate mesh, different simulations were carried out, varying the radial viscous resistance. All cases studied are shown in Table 5. Celik et al. [23] developed the grid convergence index (GCI), based on the Richardson Extrapolation, for mesh convergence analysis. This method estimates the solution by extrapolating the solutions from existing meshes and by the relative grid convergence index of the meshes produced [24]. Celik et al. [23] report that the procedure for determining the GCI starts by determining the representative mesh size, h (Equation (25)), as follows: (25) where N m is the number of mesh elements and ∆V i is the volume occupied by element i.
Using the value obtained in Equation (25), as a reference, meshes are generated, with different numbers of elements. The method determines that the ratio, r = h/h re f ined , must be greater than 1.3 for each generated mesh.
In this methodology h 1 > h 2 > h 3 , that is, h 1 will correspond to the most refined mesh and h 3 to the less refined mesh, and φ 1 , φ 2 , and φ 3 will be the respective results of a given variable analyzed. Therefore, the ratios between meshes r 21 and r 32 were defined, according to the following equations: Using Equations (28)-(30), the order of convergence, or apparent order (p), is calculated as follows: where, From Equation (29), it is observed that q(p) = 0 for r 21 = r 21 .
According to Paudel and Saenger [25], the value of the constant c that determines the convergence (Equation (33)) is given by: For c > 1 there is a monotonic divergence of the solution, 0 < c < 1 a monotonic convergence, −1 < c < 0 an oscillatory convergence, and c < −1 indicates an oscillatory divergence.
The extrapolated solutions and the approximate relative errors can be determined by Equations (34)-(36), respectively.
Therefore, based on the equations presented, the Grid Convergence Index (GCI) can be obtained using Equations (37) and (38), as follows: To obtain the radial viscous resistance (in the y-direction) corresponding to the membrane resistance used by Liao et al. [18] in their experiments, six simulations (Cases 04, 05, 06, 07, 08, and 09) were carried out, keeping constant values for feed fluxes and axial viscous resistance (1/α x ), and varying the radial viscous resistance 1/α y .
After the end of each simulation, a process parameter known as clearance was determined. This parameter represents the solute removal rate and is defined as follows: With the data obtained for each simulated case, a linear regression of the clearance parameter was performed as a function of radial viscous resistance (Equation (40)) using Excel software.
where a and b are constants to be determined. After determining the parameters a and b of Equation (40), and using the Clearance obtained experimentally by Liao et al. [18], a value of 1/α y was obtained and implemented in Ansys Fluent ® software. With these data, a new simulation was made, where a new Clearance was obtained, which was compared with the one obtained experimentally by Liao et al. [18] under the same operating conditions. Once the error of this comparison was verified, the value of 1/α y was corrected and the process was repeated until a minimum error was obtained (trial and error method).

Mesh Quality Assessment
As already mentioned, the Grid Convergence Index (GCI) method was used to evaluate the quality of the developed meshes. To perform the analyses, three computational meshes with different refinement levels were generated (M 1 , M 2 , M 3 ) under a refinement ratio of 1.5, applied to the Sizing command in the mesh generation software, Meshing by Ansys ® . This refinement ratio is per the methodology proposed by Roache [26]. The number of elements of the meshes used can be seen in Table 3. The meshes were made in a structured way, with a standardized refinement throughout the domain; details can be seen in Figures 3 and 5.
For the mesh test, the hydrodynamic variables of urea velocity and pressure were investigated using the GCI. For the analysis of urea velocity using the GCI method, three radial lines were chosen at the axial positions 20.0 mm, 101.6 mm, and 183.2 mm, as shown in Figure 6, aiming to analyze the regions close to the inlet and outlet and center of the equipment.   Table 6 presents the results obtained in the study of the GCI for the urea velocity in the y = 0.159 mm position, at the axial positions x = 20 mm, 101.6 mm, and 183.2 mm.
After analyzing Table 6, it is observed that the values of the coefficient c are in the range of 0 < c < 1, indicating monotonic convergence. Furthermore, the GCI 21 < GCI 32 indicates that the dependence of the results on the mesh elements sizes has been reduced and is approaching an independence condition. It can also be seen that the values of GCI 21 and GCI 32 are below the 10% limit, as established by Celik and Karatekin [27]. The values of the variables GCI 32 and r p GCI 21 are close enough, indicating that the extrapolated solution is close enough to the exact solution for this variable. Table 6. Parameters obtained from the study of the Grid Convergence Index for urea velocity as response variable (y = 0.159 m).

Parameter
Axial Position It is observed that, as the mesh is refined, the value of the variable of interest approaches the value of the asymptotic solution. This behavior is evidenced in the three axial positions. Meshes M 2 and M 1 present similar solutions and a slight difference is verified for mesh M 3 , with less refinement, especially in the axial position 101.6 mm.  Table 7 summarizes the obtained mean relative error of the variable compared to the obtained value using the extrapolated mesh, for the three axial positions x 1 , x 2 , and x 3 . There is a smaller mean relative error for the most refined mesh (M 1 ) compared to the extrapolated mesh, with values of 1.44%, 0.82%, and 0.2%. Similar results are observed for the mesh M 2 , with mean relative errors close to those obtained for M 1 . To ensure the quality (convergence) of the analyzed meshes and assist in decision making, in addition to urea velocity, one more hydrodynamic variable was evaluated: the pressure. Thus, the pressure was investigated using the GCI method, in the axial positions 10.0 mm, 101.6 mm, and 193.2 mm, as shown in Figure 8.  Table 8 shows the pressure values at positions x 1 , x 2 , and x 3 , as well as the parameters used in the GCI method. It is observed that the parameter c is in the range between 0 and 1, indicating monotonic convergence, with a solution within the asymptotic range, as expected for GCI 32 ∼ = r p GCI 21 [26]. Parameters GCI 21 and GCI 32 have values below 10%, satisfying the specifications for convergent solutions as proposed by Celik and Karatekin [27]. Furthermore, the GCI 21 is smaller than the GCI 32 , indicating independence of the results related to the refinement of the mesh. Analogous behavior can be observed for the other two axial positions, x 2 and x 3 .     Table 9 illustrates the mean relative error compared to that obtained with the extrapolated mesh, for the three axial positions: 10.0 mm, 101.6 mm, and 193.2 mm, for meshes M 1 , M 2 , and M 3 . Corroborating the results presented in Figure 9 it is observed that the mesh M 1 has the lowest mean relative error compared to the extrapolated mesh, with values lower than 2%, at the three positions analyzed. A similar result is observed for the mesh M 2 , which presents mean relative errors close to the values obtained for M 1 . However, the mesh M 2 presents a lower computational effort, with reduced simulation time compared to that observed for the mesh M 1 , in addition to solutions independent of the element size. It is, therefore, the best option for subsequent simulations.

Clearance
As a representative parameter of the volumetric flux of urea removed from the bloodstream, clearance was used to validate the mathematical model developed. The statistical parameters of Equation (40) achieved after adjustment to the clearance data obtained in the simulation of cases 04 to 08 are a = −10.19 mL/min and b = 559.54 mL/min. The coefficient of determination obtained after fitting was R 2 = 0.84. Figure 10 illustrates the behavior of clearance as a function of the radial viscous resistance of the membrane obtained in cases 04 to 08. It is observed that the clearance decreases with increasing viscous resistance since the permeation through the membrane included in the mathematical model is controlled by the action of the viscous resistance, responsible for regulating the transport of particles in the porous domain. Therefore, when viscous resistance is increased, there is a decrease in the flow through the membrane and, consequently, a reduction in the removal rate of urea from the bloodstream.  Table 10 summarizes the comparison between the best result obtained for clearance after the adjustment process applying the trial and error technique (Case 9) and the experimental and numerical data reported by Liao et al. [18], which were obtained under the same operating conditions. Upon analyzing Table 10, it is verified that the clearance obtained after successive simulations (Case 9) presents a value close to the experimental result reported by Liao et al. [18], with a very small error, 0.08%, much lower than the error obtained by Liao et al. [18], 6.81%. This difference is probably due to the difference in the mathematical models used in this research and the one used by Liao et al. [18]. In their research, Liao et al. [18] use a methodology based on the coupling of the domains (shell, lumen, and membrane), considering the shell domain as a porous zone, with the Darcy equations being applied to predict the flux on the shell side, the Navier-Stokes equations to simulate lumen-side flux, and the Kedem-Katchalsky equations to calculate transmembrane flux. The present work, on the other hand, uses a multiphase Eulerian-Eulerian approach for all domains, incorporating an extra resistance in the linear momentum equation to predict the flow in the porous domain, making it, therefore, a more robust and accurate model. Figure 11 shows the urea volume fraction field in the XY plane at Z = 0, for different moments of the process-500, 1000, 1500, 2000, and 2500 s-for case 09. Analyzing this figure, it is observed that the fluid takes approximately 2500 s to travel through the entire hollow fiber membrane, with higher concentrations of urea, as predicted, in the vicinity of the lumen inlet region. In the membrane, it is possible to see the axial permeation of urea due to the characteristic of the anisotropic porous medium, with a lower axial viscous resistance than the radial one (Figure 11a-c); however, this behavior is not maintained for t > 1500 s, due to the approximation of the flow to the shell inlet. Figure 12 shows the urea volume fraction field in the domains: shell, porous membrane, and lumen, in the XY plane at Z = 0, at time t = 6200 s, for case 09. It is observed a decrease in the volume fraction of urea at the end of the blood flow region and the beginning of the dialysate flow region. This behavior is associated with the countercurrent flow and the characteristics of the porous medium, which provide permeation of the solute, with higher concentrations of urea in the feed and dialysate output. Figure 13 shows the urea volume fraction profile in three axial positions-20.0 mm, 101.6 mm, and 183.2 mm-inside the porous domain, at t = 6200 s. A practically constant volume fraction profile is observed in the axial positions of 20.0 mm and 101.6 mm, and a decreasing behavior with the radial position only in the position of 183.2 mm, precisely due to the proximity to the dialysate inlet in countercurrent. Figure 14 illustrates the urea local velocity field in the domains: shell, lumen, and porous membrane, in the XY plane at Z = 0, at time 6200 s, for case 9. Upon analyzing this figure, higher local urea velocities are observed in the inlet and outlet sections of the lumen, and the exit region of the shell. Countercurrent fluxes, through a membrane with low resistance to flow and the drag force, help transport the solute, causing larger volume fractions of urea in high-velocity regions. Furthermore, variations in velocity inside the porous medium can be seen, associated with the anisotropic membrane, with axial viscous resistance (7.75 × 10 8 m −2 ) greater than the radial viscous resistance (2.15 × 10 14 m −2 ), providing greater permeation along the membrane. It is important to note that the scale of values in the legend is different. Figure 15 shows the urea velocity profile in three axial positions-20.0 mm, 101.6 mm, and 183.2 mm-inside the porous domain at t = 6200 s. Analyzing this figure, we can observe that the urea velocity decreases with the increase of radial position, due to the resistance imposed on the flow by the porous membrane. In addition, there is a similar behavior in the three positions analyzed, with a greater velocity gradient in the axial position 183.2 mm. Figure 16 represents the pressure field inside the domains of shell, lumen, and porous membrane, in the XY plane at Z = 0, at time 6200 s, for case 9. Analyzing this figure, high pressures are verified in the entrance regions of the domains, due to the mass flow in these regions, presenting a greater pressure variation in the shell region. This pressure distribution inside the equipment is associated with strict control of the dialyzer output conditions. High inlet pressures help in convective transport of solute, especially in the region close to the lumen inlet, as well as in the appearance of a possible retrofiltration in the vicinity of the lumen outlet.         Figure 18 shows the behavior of urea velocity in the form of streamlines in the domains: shell, lumen, and porous membrane, in the XY plane at Z = 0, at time 6200 s, for case 9. It is observed, by the streamlines, a behavior similar to that of Figure 14, with the greater mass flow in the region close to the lumen inlet. Furthermore, it is possible to visualize the drag of particles and solute permeation through the membrane, in addition to the presence of the retrofiltration phenomenon. Figure 19 shows the urea velocity vectors field at the fiber-hollow membrane interface in the Y direction, for case 9. It is possible to observe a higher urea permeation in the first half of the membrane, reducing with the progress of the flow due to the influence of countercurrent dialysate flow. The retrofiltration phenomenon is evident when observing the velocity field in the form of vectors; it can be seen that the fluid flows from the lumen to the shell in the left region of the domain and the opposite occurs in the right region of the domain. This behavior is expected for hemodialysis membranes, as observed by Eloot [28].

Conclusions
Based on the numerical results obtained in the simulations of the hemodialysis process via hollow fiber membrane, it can be concluded that the proposed mathematical model proved to be useful to predict the fluid behavior inside the hollow fiber membrane, allowing a better understanding of the fluid dynamics inside the equipment. The clearance parameter decreases with increasing radial viscous resistance and the analysis of the volume fraction field and streamlines of urea showed greater permeation in the vicinity of the lumen inlet region. Moreover, the highest pressures were observed at the shell and lumen inlets, with a greater variation on the shell side, due to the presence of the hollow fiber membrane separating the two domains. Further, the presence of retrofiltration phenomena was observed, especially in the vicinity of the lumen outlet region, due to the influence of countercurrent dialysate feeding.

Data Availability Statement:
The data that support the findings of this study are available upon request from the authors.