- Mathematical model development theory
The mathematical model makes following assumptions: (1) Intraparticle diffusion is considered as rate controlling step (2) The ion-exchange fiber is modeled as a slender shelled cylinder with shell thickness in the order of 10-9m, which is much smaller than the cylinder’s length, so as to justify one-dimensional diffusion taking place only in the radial direction. (3) The thickness is small compared (῀10 nm) to the radius (῀350 nm), 1-D diffusion, therefore, is modeled as diffusion within a thin flat shell. (4) The shell is uniform in size and density. (5) A partitioning exists between the outer surface of the cylinder and the bulk solution.(6) The ratios of diffusion coefficients of Na+ and Cu2+, and Ca2+ and Cu2+ within the fiber is taken to be same as the ratios of the diffusion coefficients of these ions in free form in water. It is assumed that the conditions within the fiber is similar for all the ions, and their relative mobilities are equally affected.
Starting with the generalized model considering both the Fickian diffusion term and Nernst-Plank term for both the species, and overall electro-neutrality condition, the model is rigorously derived below. The final derivation shows that the generalized model effectively reduces to a Fickian diffusion model.
The model derivation follows:
The ion-exchange taking place inside the fiber:
A, in our case, is either Na+ or Ca2+ and designating Cu2+ as B
Using the generalized model using both Fick’s law – for diffusion term, and Nernst-Plank model for influence of the electric field caused due to the difference in the mobilities of the counter-ions on the ion-exchange rate, the flux equation of species A can be written as:
The flux equation for species B
where, N = the flux of respective ions; D = the diffusivity inside the fiber; C = the concentration of respective ions inside the fiber; Z = the chemical valence of ions (2 in this case for Cu2+, and Ca2+ and 1 for Na+); F = Faraday’s constant; R = universal gas constant; T = temperature; Ø = electric potential
To maintain electro-neutrality, it is assumed that no electric current develops i.e. net flux is zero and there are no vacant ionic sites inside the fiber. Therefore, following equations satisfy the above conditions:
No current condition
No vacant sites condition
where V = Volume inside the fiber
Avogadro No. = 6.023 X 1023
Mweightrefers to the molecular weight of the respective ions
CB(t=∞) is the equilibrium concentration of the species diffusing into the fiber
For Ca2+ - Cu2+ exchange kinetics, equation 4 becomes:
Also,
For Na+ - Cu2+ exchange kinetics, equation 4 becomes:
Also,
Also, the ratios of diffusivities of Ca2+ and Cu2+ and Na+ and Cu2+ in the fiber is taken same as the ratios of their free ion diffusivities in water. The ratios are calculated from the following diffusivity value mentioned in the literature [Reference #1 of the Supporting Information]:
Diffusivity of free Ca2+ in water = 7.92 X 10-10 m2/s
Diffusivity of free Na+ in water = 1.33 X 10-9 m2/s
Diffusivity of free Cu2+ in water = 7.33 X 10-10 m2/s
Therefore, DA/DB for Ca2+ - Cu2+ exchange kinetics is taken as 1.08 and for Na+ - Cu2+ exchange kinetics as 1.8 respectively.
Solving equation 1 to 3 using equations 5 & 6 and diffusivity ratio 1.08 (for Ca2+ - Cu2+ exchange kinetics) and equations 7 & 8 and diffusivity ratio 1.8 (for Na+ - Cu2+ exchange kinetics) following expressions for Cu2+ flux into the fiber for Ca2+ - Cu2+ exchange and Na+ - Cu2+ exchange were obtained respectively:
Equation 9
Equation 10
Since, CB/CB(t=∞) varies from 0 to 1, the functions, f1 = and, f2 = vary from 0.68 to 1 and 0.65 to 1 respectively. In our case, we have taken effective diffusivity as
in the case of Ca2+ - Cu2+ exchange Equation 11
and
for Na+ - Cu2+ exchange Equation 12
Thus,
This essentially converts the Cu2+ flux Fickian with variable diffusivity. Since, in our case we were interested in establishing the value of lumped effective diffusion coefficient for the slender shell, we took it as constant, and also the variations of the functions f1 and f2 are between 0.65 to 1, and do not affect much the order of magnitude of the predicted diffusivity.
Using Fick’s law and applying mass balance, the concentration of the counter-ion satisfies the following equation:
(1)
The initial and boundary conditions of equation (1) are:
(2)
(3)
(4)
Where = concentration of solute inside the fiber; = intraparticle diffusivity; r0 and R0 are inner and outer radii respectively; K = partition coefficient; LN = the effective length of the fiber;
Parameters such as partition coefficient (K) and the value of LN were calculated as follows:
LN calculation
The fiber is modeled as a single cylinder with an effective length of LN, where N = number of fiber strands and L = length of one fiber strand. It is calculated as follows:
Where = density of Poly(acrylonitrile); = mass of fiber taken
Partition coefficient calculation
By definition of partition coefficient (K)
(7)
Using (7) and (8), partition coefficient (K) is calculated as:
The equations 1-5, with values obtained from the equations 6 and 9, are numerically solved using Crank – Nicholson method in MATLAB R2011b.
- Numerical Simulation Scheme
The governing equation on discretization following Crank-Nicolson method yields to
Where r = ; i is a subscript for spatial co-ordinate and varies in equation (10) from 2 to Nx– 1(i = 1 is left hand boundary and i = Nxis right hand boundary); k is a subscript for time co-ordinate (k = 1 is t = 0 and k = Nt is tmax).
Nxis calculated as: Nx = (xmax – xmin)/∆x + 1; Nt = tmax/∆t + 1
Spatial co-ordinates have been non-dimensionalized using x = , and D = .
At the left hand boundary, Neumann boundary condition (eq (3)) is discretized using second order Taylor series approximation
(11)
The right hand boundary, temporal Dirichlet boundary condition, is calculated at each time step by equating mass loaded on fiber and corresponding concentration loss in the remaining bulk solution, assuming the volume of bulk solution remains constant i.e.
Equation 10, where i ranges from 2 to Nx – 1, along with equations 11 and 12 leads to a matrix equation of the form AC = B;
Where A is a tri-diagonal matrix of dimension Nt X Nx – 2; C is a column vector of dimension Nx – 2 X 1. The diagonal elements of A are 2(1+r) with exception of A(1,1) being 2 + 2r/3 due to Neumann boundary condition on the left hand side, the off-diagonal elements are –r with exception of A(1,2) being -2r/3.
The scheme is simulated in MATLAB r2011b according to the algorithm depicted in Figure 1.
Error between model obtained Cbulk and experimental Cbulk is calculated as follows
% Error = (13)
This creates a column vector of errors at each time point where experimental Cbulk is measured. Using these points root mean square error (RMSE) is calculated using inbuilt MATLAB functions. The tolerance (ϵ) for RMSE error has been set at 1% error.
Until the RMSE falls within the specified tolerance (ϵ), intraparticle diffusivity is updated using equation (5).
- Desorption/Regeneration Modeling
It is assumed that as the reaction boundary progressively moves inward, the concentration of Cu2+ in the reacted shell reduces to zero, and the concentration of unreacted core remains at C0: the initial value of what it was prior to regeneration. Therefore, mass lost to the bulk, as measured at different intervals by drawing a fixed amount of aliquots, is used to determine the movement of reaction front inwards as follows:
Using mass lost versus time plot, slope ( ) is calculated for different time intervals. A time interval (2 minutes in this case) is chosen, and mass lost at this interval is calculated for the fiber. Based on these values, progress of reaction front every time step, and radius of unreacted core at end of every 2 minutes are calculated.
Say at t = t, radius has of unreacted core is shrunk to r, and at t = t+∆t, radius further shrinks by ∆th. Therefore, mass lost, dm is given by
dm = C0πLN(∆th(2r-∆th)) (14)
as 2r > ∆th, (as 2r ῀700 nm & ∆th ῀10 nm), equation 14 reduces to
dm = 2πrC0LN∆th (15)
Magnitude of the diffusive flux moving through this thickness ∆th is given by:
Therefore, the mass lost in shrinking by ∆th in ∆t time is written as:
Substituting this mass lost obtained from equation (15) as
dm = (16)
is solved for every ∆t = 2 minutes until more than 90% of the regeneration is achieved. Finally, an average value of these’s is taken and compared with the value obtained from distributed model numerical solution.
Desorption/regeneration modeling using SCM predicts intraparticle diffusivity in the order of 10-20 m2/sec. It is seen that more than 90% of regeneration is achieved in first 30 minutes (Table SI), and penetration of the reaction front is almost over within this time (9.7 nm as compared to 10 nm overall). Therefore, it is seen that overall area average diffusivity (2.5 x 10-20 m2/s) is of the same order of magnitude as in the region where around 72% of the mass leaches (3.48 x 10-20m2/s) (Table SI).
Table S1: SCM modeling result for desorption data
Time Interval (min) / Radius @start of interval (nm) / Radius @end of interval (nm) / Intraparticle Diffusivity (m2/s) / Incremental Mass Leached (%) / Cumulative Mass Leached (%) / Area AverageDiffusivity (m2/s)
0-10 / 350 / 343 / 3.48 x 10-20 / 71.6 / 71.6 / 2.5 x 10-20
10-20 / 343 / 342.1 / 3.1 x 10-22 / 9.6 / 81.2
20-30 / 342.1 / 341.1 / 3.1 x 10-22 / 9.6 / 90.8
30-40 / 341.1 / 340.7 / 1.7 x 10-23 / 2.2 / 93.0
40-50 / 340.7 / 340.5 / 1.7 x 10-23 / 2.2 / 95.2
50-60 / 340.5 / 340.3 / 1.7 x 10-23 / 2.2 / 97.4
60-70 / 340.3 / 340.23 / 6.3 x 10-25 / 0.43 / 97.83
70-80 / 340.23 / 340.18 / 6.3 x 10-25 / 0.43 / 98.26
80-90 / 340.18 / 340.14 / 6.3 x 10-25 / 0.43 / 98.69
90-100 / 340.14 / 340.1 / 6.3 x 10-25 / 0.43 / 99.12
100-110 / 340.1 / 340.05 / 6.3 x 10-25 / 0.43 / 99.55
110-120 / 340.05 / 340.01 / 6.3 x 10-25 / 0.43 / 99.98
Reference:
1.)Yuan-Hui, L., & Gregory, S. (1974). Diffusion of ions in sea water and in deep-sea sediments. Geochimica et cosmochimica acta, 38(5), 703-714.