Simulation of Ion Exchange Resin with Finite Di ﬀ erence Methods

: Ion exchange resin is used to remove potentially corrosive impurities from coolant in the ﬁrst circuit of a nuclear power plant. After one operational cycle, the used and unused resin in the mixed bed is discarded as solid waste. The aim of this work is to create a mathematical model to predict the operational cycle time of the mixed bed resin for reducing unused resin discharge. A partial di ﬀ erential equation (PDE) was set up with the conservation of matter. A ﬁnite di ﬀ erence method was used to solve the PDE. Matlab was the programming and calculating tool used in this work. The data from solution were obtained at di ﬀ erent time and space nodes. The model was then veriﬁed experimentally using di ﬀ erent ions on exchange columns. Concentrations of K + , Mn 2 + , and Cl - were calculated to verify the validation of the model by comparing it with experimental data. The calculated values showed good consistency with the experimental value.


Introduction
A chemical and volume control system is an auxiliary system in the first circuit of a nuclear power plant [1][2][3][4]. Its main function is to regulate the purity of coolant and to moderate mechanical stresses (such as pressure) on the circuit. Purification is achieved through the use of a mixture of cationic and anionic exchange resins in the fixed bed [5][6][7][8][9].
The ability to calculate breakthrough curves [10,11] for mixed bed ion exchange resins [12][13][14][15][16][17] (i.e., the relationship between running time and concentrations of pollutant ions that escape in the effluent) would be of great use in predicting the operational lifetimes of mixed-bed resins for nuclear reactors [18][19][20]. However, most of these models are limited because they are either difficult to solve or require too many dynamic parameters.
The first attempt at creating such a model was made in 1951, by J. B. Rosen of Oxford University [21]. Further work was done by Yi et al. [22], Ebrahimi. et al. [23], Xiu [24], and Wu et al. [25], who mainly based their models on equations for linear driving force and thin film diffusion [26,27] at low solute concentrations [28]. All of the above were typical models for predicting breakthrough curves, which provided some inspiration for this work.
Some of the models most commonly used today for modeling ion exchange are the Thomas (TH), Yoon-Nelson (YN), and Adams-Bohart (AB) [29,30] models. However, one common characteristic of all three models is that they produce a linear function that relates effluent concentration to running time. Further work by Kataoka [31,32] showed the influence of diffusion coefficient on computation. Lin, X. et al. [33], Aurélie et al. [34], Yang, J. et al. [35], and Murillo, R. et al. [36] also produced methods for modeling ion exchange. However, the concentration of the other position column, except outflow, could be calculated and are seldom mentioned.
In summary, all previous research has modeled ion exchange using an analytic method [37]. There are few studies that have attempted to model ion exchange using a finite difference method. Yi [22] and Jian Wu [25] used a finite difference method to solve an ion exchange model, and the difference compared to this work was that the complete derivation and discretization of the difference equation. The results of their work produced comprehensive computational data. The advantage of a finite difference method is that the results of the calculation are presented from two dimensions of time and space. A finite difference method requires much more computation than an analytical method, however the development of modern computers has solved this problem.
This work aims to use a finite difference method to produce a novel model for calculating breakthrough curves of ion exchange resins. Partial differential equations that describe changes in the concentration of effluent over time were derived using a differential element method. The partial differential equations (PDEs) were then solved using a finite difference method with controllable precision and good convergence. Calculation and data visualization were performed using Matlab. This model was then tested against experimental data and found to be in good agreement. The figures of data have intuitive visualization.

Modeling Assumptions
The system is controlled entirely by film diffusion [4]. The selectivity coefficient is constant and independent of temperature. The reaction speed is instantaneous as compared to the exchange speed. Axial diffusion can be ignored. The whole flow process is plug flow [27]. All changes to the system are isothermal and isobaric. The resin is homogeneous. There is no net electrical current.

Equations
Wolfgang [38] provided new information about the combination of the Nernst-Planck equation [39,40] and convection equation. The basic idea of this method was to discretize the ordinary equation to get the difference scheme. The detailed derivative is shown in Appendix A.

Initial Condition
On the basis of the actual situation, the concentration value at the inlet of the column is the initial concentration of solution at any time. Obviously, there is no substance when the time value is zero. Thus, the initial conditions are: The initial conditions are discretized.

Numerical Method
When calculating these models, it is necessary to compare the actual data of the test with the existing article parameters, such as the length of the resin column, velocity, and initial concentration. In general, modeling the behavior of the ion exchange resin over time begins by substituting the initial conditions into the system of PDEs and then proceeds iteratively until the concentration values at all space and time nodes are obtained.
The selection of an appropriate time step and space node is very significant for successful calculation using this system. As mentioned above, a one-dimensional convection equation is solved using a finite difference method to yield an upwind difference scheme. Selection of the time step directly affects the results.
A program was created in Matlab that used a for loop structure to calculate the concentration data at different time nodes, which proved to be an effective way to carry out the iterative calculation. Unlike an analytic solution to the system of equations, this technique produced a matrix containing concentration values at all space and time nodes for the flow of material through the column.

Experiment
The nuclear grade purification resin (ZGNR8715) was provided by the ZhengGuang company. Its physical and chemical properties are listed in Table 1. The experimental procedure is shown in Figure 1. Tests were conducted under controlled conditions of pressure and temperature. The breakthrough point was determined when the conductivity of effluent was 2 µs/cm. The water used in the experiment was ultrapure water filtered via reverse osmosis, electro-deionization, and a mixed bed ion resin purification column (separate from the resin being tested). The conductivity range of the ultrapure water was 0.06~0.16 µs/cm. A constant-speed pump was used to inject 3 × 10 −4 mol/L solution of a given ion into the top of the column. The column was filled with 100 mL of resin. The effluent samples were then collected at intervals and analyzed using high performance ion chromatography. Unlike an analytic solution to the system of equations, this technique produced a matrix containing concentration values at all space and time nodes for the flow of material through the column.

Experiment
The nuclear grade purification resin (ZGNR8715) was provided by the ZhengGuang company. Its physical and chemical properties are listed in Table 1. The experimental procedure is shown in Figure 1. Tests were conducted under controlled conditions of pressure and temperature. The breakthrough point was determined when the conductivity of effluent was 2 μs/cm. The water used in the experiment was ultrapure water filtered via reverse osmosis, electro-deionization, and a mixed bed ion resin purification column (separate from the resin being tested). The conductivity range of the ultrapure water was 0.06~0.16 μs/cm. A constant-speed pump was used to inject 3 × 10 −4 mol/L solution of a given ion into the top of the column. The column was filled with 100 mL of resin. The effluent samples were then collected at intervals and analyzed using high performance ion chromatography.

Convergence of the Model
The validation of a model not only depends on the consistency but also relies on convergence. In the circumstance that the computation is stable, values that are too big or too small have an effect on accuracy. A smaller time step can make the result more precise but can also create a risk of nonconvergence. Thus, exploring the effect of the time step value on model convergence is necessary. The effect of the interval value is shown in Figure 2.

Convergence of the Model
The validation of a model not only depends on the consistency but also relies on convergence. In the circumstance that the computation is stable, values that are too big or too small have an effect on accuracy. A smaller time step can make the result more precise but can also create a risk of non-convergence. Thus, exploring the effect of the time step value on model convergence is necessary. The effect of the interval value is shown in Figure 2. The effect of the interval value is shown in Figure 2. When the dt value is 0.1 or 0.01, there are totally different computation curves. The blue line, in Figure 2, shows a good result, but the hot line with a value of 0.01 exhibits non-convergence. By constantly adjusting the time step, the most appropriate value is determined. In a complex iterative process, the value of the time interval needs to be accurate enough to make the result adaptable. One advantage of Matlab programming is that adjusting a parameter is really easy, which means the program can be debugged when the calculation result is not converged.

Analysis of Theoretical Results
Accuracy of the model was assessed by comparing the calculated values of effluent over time to the experimentally determined values for several different ions. Visualization of data is a unique aspect in this work, and convergence of the equation is essential to a successful calculation with this model. The equation could be adjusted at the same time to achieve the purpose of getting close to the actual situation. The predicted curves were obtained by using the last space point data from the matrix and conducting iterative calculations. Experimental data for Mn 2+ were obtained using a realtime conductivity meter. It was assumed that when the value of C/C0 reached 0.05, the resin was at the breakthrough point. In general, the model closely matched the experimental data for Mn 2+ , which is of particular significance because the Mn 2+ ion is the main impurity in the coolant of a nuclear reactor. Thus, the main calculation target is based on the experimental condition of Mn 2+ . To further validate the model, calculated predictions were also compared to experimental data for Cl -, and the results are shown in Figure 3.  Figure 4. The two groups of data had a good fitting degree and a consistent trend of change for their correlation coefficients, which were 0.9939 and 0.9893, The effect of the interval value is shown in Figure 2. When the dt value is 0.1 or 0.01, there are totally different computation curves. The blue line, in Figure 2, shows a good result, but the hot line with a value of 0.01 exhibits non-convergence. By constantly adjusting the time step, the most appropriate value is determined. In a complex iterative process, the value of the time interval needs to be accurate enough to make the result adaptable. One advantage of Matlab programming is that adjusting a parameter is really easy, which means the program can be debugged when the calculation result is not converged.

Analysis of Theoretical Results
Accuracy of the model was assessed by comparing the calculated values of effluent over time to the experimentally determined values for several different ions. Visualization of data is a unique aspect in this work, and convergence of the equation is essential to a successful calculation with this model. The equation could be adjusted at the same time to achieve the purpose of getting close to the actual situation. The predicted curves were obtained by using the last space point data from the matrix and conducting iterative calculations. Experimental data for Mn 2+ were obtained using a real-time conductivity meter. It was assumed that when the value of C/C 0 reached 0.05, the resin was at the breakthrough point. In general, the model closely matched the experimental data for Mn 2+ , which is of particular significance because the Mn 2+ ion is the main impurity in the coolant of a nuclear reactor. Thus, the main calculation target is based on the experimental condition of Mn 2+ . To further validate the model, calculated predictions were also compared to experimental data for Cl -, and the results are shown in Figure 3. The effect of the interval value is shown in Figure 2. When the dt value is 0.1 or 0.01, there are totally different computation curves. The blue line, in Figure 2, shows a good result, but the hot line with a value of 0.01 exhibits non-convergence. By constantly adjusting the time step, the most appropriate value is determined. In a complex iterative process, the value of the time interval needs to be accurate enough to make the result adaptable. One advantage of Matlab programming is that adjusting a parameter is really easy, which means the program can be debugged when the calculation result is not converged.

Analysis of Theoretical Results
Accuracy of the model was assessed by comparing the calculated values of effluent over time to the experimentally determined values for several different ions. Visualization of data is a unique aspect in this work, and convergence of the equation is essential to a successful calculation with this model. The equation could be adjusted at the same time to achieve the purpose of getting close to the actual situation. The predicted curves were obtained by using the last space point data from the matrix and conducting iterative calculations. Experimental data for Mn 2+ were obtained using a realtime conductivity meter. It was assumed that when the value of C/C0 reached 0.05, the resin was at the breakthrough point. In general, the model closely matched the experimental data for Mn 2+ , which is of particular significance because the Mn 2+ ion is the main impurity in the coolant of a nuclear reactor. Thus, the main calculation target is based on the experimental condition of Mn 2+ . To further validate the model, calculated predictions were also compared to experimental data for Cl -, and the results are shown in Figure 3.  Figure 4. The two groups of data had a good fitting degree and a consistent trend of change for their correlation coefficients, which were 0.9939 and 0.9893,  Figure 4. The two groups of data had a good fitting degree and a consistent trend of change for their correlation coefficients, which were 0.9939 and 0.9893, respectively.
Thus, it can be concluded that the model verified the test data. Breakthrough points were obtained, as shown in Figure 4. As a practical matter, the breakthrough point can be used to determine when the resin used in the first circuit of a nuclear reactor should be discarded. respectively. Thus, it can be concluded that the model verified the test data. Breakthrough points were obtained, as shown in Figure 4. As a practical matter, the breakthrough point can be used to determine when the resin used in the first circuit of a nuclear reactor should be discarded.  Table 2. The comparative results between the experimental value and the calculated value meet our modeling requirements. When calculating different ions, the modification coefficient of the model is very close, indicating that the model is stable.
The advantage of a finite element method is the integrity of the solution, and thus the concentration changing processes are shown in Figures 5 and 6. Previous articles just exhibited the concentration changing curve in the effluent. Figure 5 shows the different changing processes from inflow to outflow. Putting every concentration changing curve in the same picture can directly show the difference in different positions.  Figure 5 shows the concentration curves calculated every 100 space intervals in the column. Starting from the entrance, 11 curves of concentration over time were obtained by calculating every 100 spatial nodes. It is clear that the concentration changes with time at different spatial points. In particular, from the inlet to outlet of the column, the slope decreases because it is the same as the actual situation. Obviously, the ion exchange process, indeed, is fiercer because the local capacity of  Table 2. The comparative results between the experimental value and the calculated value meet our modeling requirements. When calculating different ions, the modification coefficient of the model is very close, indicating that the model is stable.
The advantage of a finite element method is the integrity of the solution, and thus the concentration changing processes are shown in Figures 5 and 6. Previous articles just exhibited the concentration changing curve in the effluent. Figure 5 shows the different changing processes from inflow to outflow. Putting every concentration changing curve in the same picture can directly show the difference in different positions. respectively. Thus, it can be concluded that the model verified the test data. Breakthrough points were obtained, as shown in Figure 4. As a practical matter, the breakthrough point can be used to determine when the resin used in the first circuit of a nuclear reactor should be discarded.  Table 2. The comparative results between the experimental value and the calculated value meet our modeling requirements. When calculating different ions, the modification coefficient of the model is very close, indicating that the model is stable.
The advantage of a finite element method is the integrity of the solution, and thus the concentration changing processes are shown in Figures 5 and 6. Previous articles just exhibited the concentration changing curve in the effluent. Figure 5 shows the different changing processes from inflow to outflow. Putting every concentration changing curve in the same picture can directly show the difference in different positions.  Figure 5 shows the concentration curves calculated every 100 space intervals in the column. Starting from the entrance, 11 curves of concentration over time were obtained by calculating every 100 spatial nodes. It is clear that the concentration changes with time at different spatial points. In particular, from the inlet to outlet of the column, the slope decreases because it is the same as the actual situation. Obviously, the ion exchange process, indeed, is fiercer because the local capacity of  Figure 5 shows the concentration curves calculated every 100 space intervals in the column. Starting from the entrance, 11 curves of concentration over time were obtained by calculating every 100 spatial nodes. It is clear that the concentration changes with time at different spatial points. In particular, from the inlet to outlet of the column, the slope decreases because it is the same as the actual situation. Obviously, the ion exchange process, indeed, is fiercer because the local capacity of resin is lower in the lower part of column. In the essence of the algorithm, the iterative calculation gives the equation the ability to exhibit the changing process of the adsorption capacity.
To observe the concentration changing process of the column more intuitively, the data of the matrix is transformed to the colorful diagram shown in Figure 6. Using imagesc function in Matlab, the concentration data is expressed using different colors. With the colorbar function, the colorbar nearing the colormap figure shows the value changing process. Thus, there are six calculation results for the column transformation in every 2000 s. resin is lower in the lower part of column. In the essence of the algorithm, the iterative calculation gives the equation the ability to exhibit the changing process of the adsorption capacity.
To observe the concentration changing process of the column more intuitively, the data of the matrix is transformed to the colorful diagram shown in Figure 6. Using imagesc function in Matlab, the concentration data is expressed using different colors. With the colorbar function, the colorbar nearing the colormap figure shows the value changing process. Thus, there are six calculation results for the column transformation in every 2000 s. As seen in Figure 6, through the built-in imagesc function in Matlab, the variation chart of concentration in the column is visualized with a very intuitive interval of 2000 s. Setting equal time interval points effectively shows the variation trend of the concentration in the column at different positions. In the Figure 6, the yellow part represents the "failure layer", the blue part represents the "inactive layer", and between them is the "working layer". This mode of visualizing the data very intuitively shows that as ion exchange continues over time, the "working layer" moves down the column. Moreover, one can see an overall trend that effluent concentration decreases over time.

Conclusions
Ion exchange resins are crucial to the functioning and safety of modern water-cooled nuclear reactors because they remove metal ions that would otherwise corrode a reactor's coolant system. However, to date there has not been a satisfactory mathematical model to predict the optimal operational lifetime of a resin bed. This work modeled resin ion exchange behavior using a system of PDEs based on the conservation of matter and then solved it using the method of finite differences. Using Matlab, the behavior of a resin bed can then be calculated via an iterative approach not merely over time, as in past studies, but over the spatial coordinate of the position of resin beads in the column. Theoretical results of this model proved to be in good agreement with experimental results using real ion exchange columns and Mn 2+ , which is the main contaminant found in reactor coolant water. Furthermore, simulation of the dynamic process is obtained, which clearly exhibits a concentration changing trend. The results show that this model has great potential to be used in nuclear engineering applications because a more efficient use of resin will not only cut costs but also reduce the quantity of toxic waste produced.  As seen in Figure 6, through the built-in imagesc function in Matlab, the variation chart of concentration in the column is visualized with a very intuitive interval of 2000 s. Setting equal time interval points effectively shows the variation trend of the concentration in the column at different positions. In the Figure 6, the yellow part represents the "failure layer", the blue part represents the "inactive layer", and between them is the "working layer". This mode of visualizing the data very intuitively shows that as ion exchange continues over time, the "working layer" moves down the column. Moreover, one can see an overall trend that effluent concentration decreases over time.

Conclusions
Ion exchange resins are crucial to the functioning and safety of modern water-cooled nuclear reactors because they remove metal ions that would otherwise corrode a reactor's coolant system. However, to date there has not been a satisfactory mathematical model to predict the optimal operational lifetime of a resin bed. This work modeled resin ion exchange behavior using a system of PDEs based on the conservation of matter and then solved it using the method of finite differences. Using Matlab, the behavior of a resin bed can then be calculated via an iterative approach not merely over time, as in past studies, but over the spatial coordinate of the position of resin beads in the column. Theoretical results of this model proved to be in good agreement with experimental results using real ion exchange columns and Mn 2+ , which is the main contaminant found in reactor coolant water. Furthermore, simulation of the dynamic process is obtained, which clearly exhibits a concentration changing trend. The results show that this model has great potential to be used in nuclear engineering applications because a more efficient use of resin will not only cut costs but also reduce the quantity of toxic waste produced.
Author Contributions: Y.Z. for writing-original draft: and software; B.L. for data curation; R.P. for review&editing; Y.L. and P.Y. for supervision.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
Based on the conservation of matter, an infinitesimal distance from the ion exchange column was studied. Flux through the column (mol/m 2 s) was defined as consisting of the mainstream and back mixing: where J is the flux of the column, u is the velocity of flow, C is the concentration at height H, and E is the back-mixing coefficient. According to the Taylor expansion: Equation (A1) with the Taylor expansion deformation can be written as: Because of the conservation of matter, one can assume that the inflow is equal to the effluent plus the accumulation in the solid and liquid phases. This gives: where, εS(J H − J H+dH ) is the inflow minus outflow in the interface εS ∂C ∂H dH is the change of ions in the microelement, and (1 − ε)S ∂q ∂t dH is the ion change in the resin.
Combining Equations (1), (3), and (4) yields: Most prior studies on this subject were based on the one-dimensional convection equation, and the various groups made different choices when deciding the magnitude of the change in the rate of adsorption. Most of them also chose to omit the back-mixing term to simplify the calculation. This assumption is indeed suitable for lab scale experiments. Therefore, the righthand side of Equation (A4) is assumed to be negligibly small: According to previous literature [18], thin-film diffusion is governed by the concentration difference equation: where d p is the diameter of the resin particles, of which the average value is given in Table 1. By combining the above two equations, the one-dimensional convection equation in this special format can be solved using the finite difference method, and the concentration values at different time and space nodes can be calculated.
The following is the derivation process of J i : Because of the assumption that there is no net electric current, the following can be written: By both sides of Equation (A9) by z i and taking the sum, the following can be derived: The total equivalent concentration c g is introduced to give: The mean value for the co-ion valence is introduced to yield: Equations (A11) and (A12) are substituted into Equation (A10) which results in the following equation: Equation (A13) is substituted into Equation (A8), and the following is obtained: where n i = − z i z Y . ξ is derived from both sides of Equation (A14) as follows: Both sides of Equation (A15) are summed to give: According to Equation (A11): This is substituted into Equation (A16) to give: According to the chain rule, the following can be written: Equations (A17) and (A18) are substituted into Equation (A15), and the following is obtained: Equation (A19) is a second-order homogeneous linear differential equation (Euler equation). According to the relationship between c g and t: c g = e t t = ln c g The following is determined by boundary conditions: For a i , the following can be written for b i : The derivative of Equation (A22) gives: The above equation is substituted into Equation (A14) to give: with P= n i : Both sides are multiplied by z i , and the equation is combined with Equation (A10) to give: With the definitions of a i and b i , the following can be written: After algebraic transformation, the following is obtained: Derivatives of Equation (A26) on both sides of ξ give: Thus, the following is obtained: Kataoka [32] mentioned the following: Both sides of Equation (A24) are multiplied by z i with the assumption that there is no net electric current and that dc g dξ is not zero, and thus the following is obtained: Equation (A31) is valid at any point in the liquid film. The value of c g is independent of the position, and thus the expressions in parentheses are all zero. n i=1 (n i + 1)D i a i = 0 (A32) For the total equivalent concentration, the following can be written: When this is combined with the definition of b i , the following is obtained: With Equations (A22) and (A26), we obtain the final target function: This PDE describing the change of concentration in the column, as obtained above using material balance and Taylor's expansion formula, is the theoretical basis and main formula in this work.
To solve this PDE using the finite difference method, upwind difference scheme of the PDE follows: This upwind difference scheme constitutes one solution to this set of PDE. It has two forms, one contains the forward difference, and the other contains the backward difference. Obviously, U L is greater than zero, and thus it is shown as Equation (38).