In [7]:
import sys
sys.path.append('../')
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
from porousmedialab.column import Column
In [8]:
tend =3
dt = 0.0025
dx = 0.2
L = 10
phi=0.99
ftc = Column(L, dx, tend, dt, w=5)
In [10]:
C_init=0.0000001/phi/3
ftc.add_species(theta=phi, name='H2CO3', D=50, init_conc=C_init, bc_top_value=0, bc_top_type='flux', bc_bot_value=0, bc_bot_type='flux')
ftc.add_species(theta=phi, name='HCO3', D=50, init_conc=C_init, bc_top_value=1e-3, bc_top_type='dirichlet', bc_bot_value=0, bc_bot_type='flux')
ftc.add_species(theta=phi, name='CO3', D=50, init_conc=C_init, bc_top_value=0, bc_top_type='flux', bc_bot_value=0, bc_bot_type='flux')
ftc.add_acid(species=['H2CO3', 'HCO3', 'CO3'], pKa=[3.6, 10.32])
ftc.add_species(theta=phi, name='H2SO4', D=50, init_conc=C_init, bc_top_value=0, bc_top_type='flux', bc_bot_value=0, bc_bot_type='flux')
ftc.add_species(theta=phi, name='HSO4', D=50, init_conc=C_init, bc_top_value=0, bc_top_type='flux', bc_bot_value=0, bc_bot_type='flux')
ftc.add_species(theta=phi, name='SO4', D=50, init_conc=C_init, bc_top_value=1e-6, bc_top_type='dirichlet', bc_bot_value=0, bc_bot_type='flux')
ftc.add_acid(species=['H2SO4', 'HSO4', 'SO4'], pKa=[-10, 1.99])
ftc.add_species(theta=phi, name='Ca', D=50, init_conc=1e-9, bc_top_value=1e-2, bc_top_type='dirichlet', bc_bot_value=0, bc_bot_type='flux')
ftc.add_species(theta=1-phi, name='CaCO3', D=5, init_conc=0., bc_top_value=0.0, bc_top_type='flux', bc_bot_value=0.00, bc_bot_type='flux', w=0.)
ftc.add_species(theta=1-phi, name='CaSO4', D=5, init_conc=0., bc_top_value=0.0, bc_top_type='flux', bc_bot_value=0.00, bc_bot_type='flux', w=0.)
ftc.add_ion(name='Ca', charge=2)
In [11]:
ftc.constants['Ks_CaCO3'] = 3.3e-9
ftc.constants['Ks_CaSO4'] = 10**-4.58
ftc.constants['k_pre'] = 1e-4
ftc.constants['k_dis'] = 1e-3
ftc.constants['CF'] = (1-phi)/phi ### conversion factor
In [12]:
ftc.rates['R_pre_CaCO3'] = 'k_pre * (Ca*CO3/Ks_CaCO3-1)'
ftc.rates['R_pre_CaSO4'] = 'k_pre * (Ca*SO4/Ks_CaSO4-1)'
ftc.rates['R_dis_CaCO3'] = 'k_dis * CaCO3 * (1 - Ca*CO3/Ks_CaCO3)'
ftc.rates['R_dis_CaSO4'] = 'k_dis * CaSO4 * (1 - Ca*SO4/Ks_CaSO4)'
In [ ]:
ftc.dcdt['CaCO3'] = 'R_pre_CaCO3/CF - R_dis_CaCO3'
ftc.dcdt['Ca'] = '-R_pre_CaCO3 + R_dis_CaCO3*CF - R_pre_CaSO4 + R_dis_CaSO4*CF'
ftc.dcdt['CO3'] = '-R_pre_CaCO3 + R_dis_CaCO3*CF'
ftc.dcdt['CaSO4'] = 'R_pre_CaSO4/CF - R_dis_CaSO4'
ftc.dcdt['SO4'] = '-R_pre_CaSO4 + R_dis_CaSO4*CF'
In [ ]:
ftc.solve()
In [ ]:
ftc.plot_contourplots()
In [ ]:
ftc.plot_saturation_index('Ca','CO3','Ks_CaCO3')
In [ ]:
ftc.plot_saturation_index('Ca','SO4','Ks_CaSO4')
In [ ]: