CSTR model (Solving an ODE in PHREEQC)

Conceptual Models > Design of conceptual models

(1/3) > >>

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.

Navigation

[0] Message Index

[#] Next page

Go to full version