-
Notifications
You must be signed in to change notification settings - Fork 22
Description
Hi, professor,
when I learn this case (Kinetic Dissolution of Quartz), I don't know how to get the value of quartz_dissolved.
When I add the command "print(quartz_dissolved)" in the "ratefun()" function, I can get a lot of value (just one time_step). However, this value is not the same with the value I got from the customer software.
I hope to get your help.
best wish!!
%pylab inline
import phreeqpython
from scipy.integrate import odeint
pp = phreeqpython.PhreeqPython('phreeqc.dat')
def ratefun(sol, quartz_dissolved, m0, A0, V):
m = m0-quartz_dissolved
rate = (A0/V)(m/m0)**0.6710**-13.7*(1-sol.sr("Quartz"))
print(quartz_dissolved) # I try this method, but it is wrong.
return rate * 1e3
solution1 = pp.add_solution({})
t = np.array([])
y = []
year = 365243600
for time, sol in solution1.kinetics('SiO2',
rate_function=ratefun,
time=np.linspace(0,5*year, 2),
m0=158.5,
args=(23.13,0.16)):
t = np.append(t, time)
y.append(sol.total_element('Si', units='mmol'))