CSTR model (Solving an ODE in PHREEQC)

[nickvz]

Hello

I have been struggling to understand how PHREEQC solves differential equations in the RATES block. I am trying to replicate a code that was published in the article by Tiruta-Barna ( Using PHREEQC for modelling and simulation of dynamic leaching tests and scenarions). The author models a stirred tank reactor by solving the following equation in the RATES block:

dC'/dt = (Q/V_leachate*(c'_in- c'))

This equation describes the change in concentration of over time for a CSTR. Initially the reactor contains 0.001 moles of reactant, in this case NaCl. Pure water enters the reactor at flow rate Q, and the solution containing NaCl exist the reactor at flow rate Q. The analytical solution to the above mentioned equation is the following:

c'= c'_0 exp(-tQ/V_leachate)

The PHREEQC code is as follows:

RATES
   FLOW_Na
-start
10     total_element = TOT("Na")
20     IF (total_element <= 0) THEN GOTO 60
30     rates = parm(1) /parm(2) * (parm(3) -total_element)
40     moles = -rate*TIME
50     IF (total_element+moles < 0) THEN moles = -total_element
60     SAVE moles
-end

SOLUTION 1
   

KINETICS 1
   FLOW_Na
-formula      Na+      1
-m      1
-m0      1
-parms      0.001      1         0 # parm(1) = Q, parm(2) = V_leachate, Parm(3) = c"_in
-tol      1e-008
-steps          3600 in 10 steps
-step_divide   1
-runge_kutta    3

USER_GRAPH 1    Compare PHREEQC and analytical solution
-axis_titles "time" " c' mol/l"
-axis_scale x_axis 0   4000
-axis_scale y_axis 0   2e-03
-initial_solutions true
-start
10 plot_xy total_time   tot("Na"), symbol_size = 0, color = Red
-end
END

The amount of Na+ that is in solution stay unchanged throughout the simulation. Am I missing something here? Please advise.

Kind regards
Nicolaus
       

[dlparkhurst]

Try the following. The changes I made were (1) changed "rates" to "rate" in line 30, (2) removed negative in line 40 (because parm(3) - total_element) was already negative, (3) removed the IF statements (optional), the IF statements can add discontinuities in the rate expression, (4)  added 1 mmol NaCl to solution, (5) made the formula NaCl in kinetics, otherwise it is equivalent to removing sodium metal, which causes redox reactions, the plus sign is ignored, and (6) added a comma in line 10 of USER_GRAPH.


RATES
FLOW_Na
-start
10    total_element = TOT("Na")
#20    IF (total_element <= 0) THEN GOTO 60
30    rate = parm(1) /parm(2) * (parm(3) -total_element)
40    moles = rate*TIME
#50    IF (total_element+moles < 0) THEN moles = -total_element
60    SAVE moles
-end

SOLUTION 1
   Na   1
   Cl   1

KINETICS 1
FLOW_Na
-formula      NaCl    1
-m    1
-m0    1
-parms    0.001    1        0 # parm(1) = Q, parm(2) = V_leachate, Parm(3) = c"_in
-tol    1e-008
-steps          3600 in 10 steps
-step_divide   1
-runge_kutta 3

USER_GRAPH 1    Compare PHREEQC and analytical solution
-axis_titles "time" " c' mol/l"
-axis_scale x_axis 0   4000
-axis_scale y_axis 0   2e-03
-initial_solutions true
-start
10 plot_xy total_time,  tot("Na"), symbol_size = 0, color = Red
-end
END

[nickvz]

Hello

Thank you, I did not notice all the simple mistakes. It must be the heat. I am in Mozambique at the moment and is hot as hell. Thank you very much for the help. I appreciate the corrections.

Kind regards
Nicolaus

[nickvz]

Hello

I just want to make sure about the addition of NaCl to the SOLUTION block. I was under the impression if you define initial moles in the kinetic block it PHREEQC reads it as the amount of moles present in the system. Is this not the case?

Kind regards
Nicolaus

[dlparkhurst]

Yes, you do define the number of moles of a kinetic reaction in the KINETICS data block, particularly if it is a certain amount of a solid that could react completely. However, your rate expression requires Na in solution to have a nonzero rate, and, for your case, the amount of the kinetic reactant is irrelevant. It is the amount of Na in solution that drives the reaction.

[nickvz]

Hello David

Can you please advise on the following:

The link below is to an example from Tony's website. It is to solve a general ODE. I am just
having some trouble to plot the variable.

http://www.hydrochemistry.eu/exmpls/RK.html

Please correct me if I am not understanding this correctly. To solve any ODE using the RATES block in PHREEQC, the following applies:

dy/dt = ((1-C)/y)*D, or in the case of the example  dy/dx = x * exp(-x)

The constants are defined in either the KINETICS block or the RATES block. The expression is integrated by multiplying the expression by TIME.

How do you plot this variable as a function of time etc. to calculate for instance the conversion as a function of time for a plug flow reactor.

Kind regards
Nicolaus

[dlparkhurst]

I think you want the functions KIN and KIN_DELTA, which give the total amount of a kinetic reactant and the change in the amount of a kinetic reactant during the last rate integration.

[nickvz]

Hello David’

I am trying to model the shrinking core model equation

d(r) / d(t) = -(((1-por)*D2)/(eps*rho_p*f_p*(r^2)))*((R*r)/(R-r))*(C_O2/H)

I have solved the ODE using the software Polymath. I used the following input for PHREEQC:

RATES
SCM
10 por      =   0.3   # Porosity
20 D2      =   2.6e-8   # Ash layer diffusion coefficient
30 eps      =   1.75   # ratio of O2 consumed to moles pyrite oxidized
40 rho_p   =   4.8   # g/cm3
50 f_p      =   0.03   # Weight fraction of pyrite in waste rock
60 R       =   0.01   # Initial core radius (cm)

70 dt = time
80 r = 0.0099
90 r = (((1-por)*D2)/(eps * rho_p * (r^2)))*((R*r)/(R-r))*(TOT("O2"))
91 IF r= 0 then goto 100
100 Save r*dt
-end

END

SOLUTION 1
Na 1
Cl 1
EQUILIBRIUM_PHASES 1
O2(g)               -0.67

KINETICS
SCM
-formula      NaCl 1
-m 1
-m0 1

-time_step 15000 in 50
-step_divide 1
INCREMENTAL_REACTIONS true

USER_GRAPH 1    Compare PHREEQC and analytical solution
-axis_titles "time" " core radius"
-axis_scale x_axis 0 15000
-axis_scale y_axis 0 0.01
-headings PHREEQC analytical
-initial_solutions true
-start

10 plot_xy total_time ,  r, symbol_size = 0, color = Red
-end
END

The following error message is shown in the output:

WARNING: Zero divide in BASIC line
 90  r = (((1-por)*D2)/(eps * rho_p * (r^2)))*((R*r)/(R-r))*(TOT("O2")).

I am guessing r is zero in the expression. I tried to correct this by adding line 80 into the RATEs block as an initial condition, but the problem is persisting so I am assuming that PHREEQC does not read this as an initial condition.

Can you please assist me with the input file and comment as to how initial conditions is defined for r(0) where r(0) = R.

Kind regards
Nicolaus

[dlparkhurst]

Variables are case independent, so R - r is zero. Also, TOT("O(0)") is the total moles of oxygen valence zero, which is one half MOL("O2"). TOT("O2") is zero.

You can use PRINT Basic statements to debug your RATES equations. The PRINT statements write to the output file. It is not pretty, but you can determine what is going wrong with your RATES equations.

[nickvz]

Hello David

I have a question about the KIN function. Does this function recall the amount of reactant (in the solid phase) that is present in the system. I need to define a value for r (uncreated core radius) of pyrite before I can solve the differential equation for the rate of decrease of r, see code below:

RATES
Pyrite
-start
10    por      =   0.3   # Porosity
20    D2      =   3e-8   # Ash layer diffusion coefficient
30    eps      =   1.75   # ratio of O2 consumed to moles pyrite oxidized
40    rho_p   =   30   # kg/m3
50    f_p      =   0.03   # Weight fraction of pyrite in waste rock
60    r_0       =   10   # Initial core radius,m
70   PI      =   3.14159
80   M_pyite   =   119.96   # g/mol
90   u0      =   0.265
100   H      =   26.3



110   r = (((3*KIN("Pyrite")*M_pyrite)/(4*PI*rho_p))^(1/3))


120   r = (-((1-por)*D2)/(eps * rho_p * (r^2)))*((r_0*r)/(r_0-r))*(u0/H)      
130   moles = (((4/3)*PI*(r^3)*rho_p)/M_pyrite)

140   save moles
-end
END


SOLUTION 1

EQUILIBRIUM_PHASES 1
O2(g)               -0.67

KINETICS
Pyrite
-formula FeS2   1
-m    0.083
-m0    0.083
-time_step 15000 in 100 steps
-step_divide 1
INCREMENTAL_REACTIONS true

USER_GRAPH 1    Compare PHREEQC and analytical solution
-axis_titles "time" " core radius"
-axis_scale x_axis 0 15000
-axis_scale y_axis 0 0.01
-headings PHREEQC SCM
-initial_solutions true
-start

10 plot_xy total_time , (((3*KIN("Pyrite")*M_pyrite)/(4*PI*rho_p))^(1/3)), symbol_size = 0, color = Red
-end
END


I have tried to define the current value for r by using the KIN function. This seems not to be working as I still have errors with division by 0.

Kind regards
Nicolaus

[dlparkhurst]

You misspelled M_pyrite in line 80.

Also, all Basic variables are local to a Basic script. Defining M_pyrite in a RATES definition does not carry over to a USER_GRAPH definition. You must either (1) define all of the needed variables in USER_GRAPH (M_pyrite, PI, rho_p), or (2) you can use PUT and GET to save data in RATES and use it in USER_GRAPH. I think case (1) would be better because RATES is evaluated many times per time step, whereas USER_GRAPH is evaluated only at the end of the time step, and values saved in RATES may actually not be the values at the end of the time step (there is an interpolation step after RATES has been evaluated).

[nickvz]

Hello
It seems like I have managed to get the model to work. I am just struggling with the diffusion coefficient. When I increase it it slows down the reaction to a near zero rate. A slower rate is expected but not this slow. The initial radius seems right at 0.00143 meter. Why does it take PHREEQC 195 seconds before it reaches this value?
Please see below is the input:

PRINT
-reset false
-kinetics
-totals


RATES
Pyrite
-start
#--------------------------------------------------------
# Define input parameters
#--------------------------------------------------------
10   por      =   0.4      # Porosity
20   D2      =   2.6e-3      # Ash layer diffusion coefficient
30   eps      =   1.75      # ratio of O2 consumed to moles pyrite oxidized
40   rho_p      =   4800      # kg/m3
50   rho_bulk   =   1500      # kg/m3
51   v_t      =   1      # bulk_liter
52   f_a      =   2/3   
60   PI      =   3.14159
70   M_pyrite   =   (119.96/1000)   # kg/mol
#--------------------------------------------------------
# Ensure r is not equel to 0 at step 1
#--------------------------------------------------------
100   IF (m < 0 ) then goto 250   
110   ratio = 1
120   if m0 > 0 then ratio = m/m0 # initial moles of the kinetic reactant
130   if ratio = 0 then ratio = 1

#--------------------------------------------------------
# Calculate volume fraction
#--------------------------------------------------------
140   mass_pyrite   = m*M_pyrite         # kg pyrite
150   v_p      = (mass_pyrite/rho_p)      # V_pyrite (l) = 0.009375 l at step 1
160   vol_frac   = v_p/v_t   # = 0.009375 l at step 1

#--------------------------------------------------------
# Calculate the initial radius of the unreacted core =r_0
# Calculate the initial radius of the particle =r_c
#--------------------------------------------------------
170   v_par   =   v_p/633.145               # Volume of individuel pyrite particles
180   r_0   =   (((3*v_par*0.857)/(4*PI))^(1/3))
190   r_c    =    r_0/(1-0.05)               # r_c is the initial core radius plus 5% ash layer
200   r = r_0*(m/m0)                     
210   if r <= 0 then goto 250

#--------------------------------------------------------
# Calculate the reaction surface in the bulk volume
#--------------------------------------------------------
220   react_a      = ((3*(10^-3))*(vol_frac/r))   
#--------------------------------------------------------
# Calculate the rate of pyrite oxidation (mol/liter_bulk.s)
#--------------------------------------------------------
230   rate = (10^3)*react_a*(D2/3)*(r_c/((r_c-r)*r))*(mol("O2"))   # mol/l_bulk.s
240   moles = rate*time
250   save moles
-end
END
#--------------------------------------------------------
# Define leachate volume
#--------------------------------------------------------

SOLUTION 1
-water      0.3
#--------------------------------------------------------
# Define EQUILIBRIUM_PHASES 1
#--------------------------------------------------------

EQUILIBRIUM_PHASES 1
CO2(g)               -2
O2(g)               -0.67

#--------------------------------------------------------
# Define KINETICS 1 for Pyrite
#--------------------------------------------------------
KINETICS 1
Pyrite
-formula FeS2   1
-m    0.375
-m0    0.375
-time_step 19500 in 100 steps
-step_divide 1
INCREMENTAL_REACTIONS true

USER_GRAPH 1    Compare PHREEQC and analytical solution
-axis_titles "time" " core radius"
-axis_scale x_axis 0 19500
-axis_scale y_axis 0 0.00146
-headings PHREEQC SCM
-initial_solutions true
-start
#--------------------------------------------------------
# Define parameters
#--------------------------------------------------------
1   rho_p         =   4800   
2   M_pyrite      =   (119.96/1000)
3   PI         =   3.14159
10   t_mol         = KIN("pyrite")
#--------------------------------------------------------
# Calculate the un-reacted core radius = r_unreacted
#--------------------------------------------------------
20   vol_pyrite      = (t_mol*M_pyrite)/rho_p
30   vol_par         = vol_pyrite/633.145
40   r_unreacted      = (((3*vol_par*0.857)/(4*PI))^(1/3))   
#--------------------------------------------------------
# Plot unreacted core radius
#--------------------------------------------------------
50    plot_xy total_time , r_unreacted, symbol_size = 0, color = Red
# 80 plot_xy total_time , kin("pyrite"), symbol_size = 0, color = Red
-end
END

[dlparkhurst]

195 s is the first time step. If you increase the number of steps or decrease the length of time, you will see a smoother transition.

I think you will have to investigate the rate equations yourself. You can use PRINT statements to look at variables to try to understand why the rate is decreasing as it is.

[nickvz]

Hello

Thank you very much. I appreciate the help.

kind regards

[nickvz]

Hello

I found the problem. Was an input error. Caused the rate to be 3 orders smaller than it is suppose to be. Thank you for all the help.

Kind regards
Nicolaus