kinetic dissolution of multicomponent glass

Conceptual Models > Kinetics and rate controlling factors

(1/1)

OLIVIA_RH:
Hello, I have written a script to model the dissolution of a glass as an alkaline solution flows through it. In order to do this, I added a new phase to the LLNL database called MultiGlass. I then input the rate equation for dissolution of this phase and input the parms into the kinetic block, in addition to initial moles of the glass (which has a GFW of ~85 g/mol).

Currently the script is failing to converge but i cannot trace the route cause of this. In the manual, i see that "no convergence" may reflect very small concentrations in the system. I have iterated through different values for m0 and parm(10) but this has not resolved the issue. Although i use these parameter values in other softwares and they return a result, I think there must be something wrong with the way i have computed the rate equation in this example.

Any help or direction on resources to help fix this would be greatly appreciated.

Thank you,
Olivia
--- Code: ---DATABASE C:\iCP\ORH Example - Silica Dissolution\phreeqc\LLNL_mod.dat

SOLUTION 0 #boundary water
temp 100
pH 11.0
units mol/kgw
density 1
Si 1e-6
Cl 1e-6
Na 0.1 charge
-water 1 # kg

SOLUTION 1-40
temp 100
pH 7
pe 4
redox pe
units mol/kgw
density 1
Na 1e-6
Cl 1e-6
Si 1e-6
-water 1 # kg

#EQUILIBRIUM_PHASES 1-40
# MultiGlass 0.0 3.5 dissolve_only

KINETICS 1-40
MultiGlass
#-formula
#-m
-m0 1
-parms 4.03e-4 0.105 45600 80000 0.30 -0.4 8.314 1 1e3 0.085

RATES
MultiGlass
# parms(1): Aa = 4.03e-4 mol/m2/s
# parms(2): Ac = 0.105 mol/m2/s
# parms(3): Eaa = 45600 J/mol
# parms(4): Eac = 80000 J/mol
# parms(5): na = 0.30
# parms(6): nc = -0.4
# parms(7): R = 8.314 J/mol/K
# parms(8): sigma = 1
# parms(9): spec As = m2/kg
# parms(10): mass of MultiGlass = kg - for 1 mol = 85 g
-start
1 Aa = parm(1)
2 Ac = parm(2)
3 Eaa = parm(3)
4 Eac = parm(4)
4 na = parm(5)
6 nc = parm(6)
7 Rg = parm(7)
8 sigma = parm(8)

9 T = TK # convert ?C to K
10 Hydrogen = ACT("H+")

# pass specific surface in m?/kg:
11 sa = parm(9)

# Compute the saturation ratio for MultiGlass:
13 qok = SR("MultiGlass")

# Calculate the individual rate contributions:
14 k_a = Aa * (Hydrogen^na) * EXP(-Eaa/(Rg*T))
15 k_c = Ac * (Hydrogen^nc) * EXP(-Eac/(Rg*T))
16 rateplus = k_a + k_c

17 mass_rock = parm(10)
18 rate = rateplus * (1 - qok^(1/sigma)) * sa * mass_rock * 0.70 * GFW("SiO2") * 1e-3
#19 save rate
20 save rate * TIME
-end

TRANSPORT
-cells 40
-length 0.002
-shifts 120
-time_step 720.0
-flow_direction forward
-boundary_cond flux flux
-diffc 0.0e-9
-dispersivity 0.002
-correct_disp true
-punch 40
END

--- End code ---
and the new phase MultiGlass was added to LLNL.dat:

--- Code: ---MultiGlass
SiO2Al0.1326504481Fe0.0022Fe0.0257Mn0.0037772087Ca0.0027055058Na0.1864276569K0.0765685019(PO4)0.0001(OH)0.7316129321 + 0.7316129321H+ = SiO2(aq) + 0.1326504481Al+3 + 0.0022Fe+3 + 0.0257Fe+2 + 0.0037772087Mn+2 + 0.0027055058Ca+2 + 0.1864276569Na+ + 0.0765685019K+ + 0.0001PO4-3 + 0.7316129321H2O
log_k 6.6731
-delta_H -2.297 kJ/mol # Calculated enthalpy of reaction #from where
# Enthalpy of formation:
-analytic -623.7513 -0.0768 41744.3699 219.6997 -2695700.6140
# -Range: 0-300

--- End code ---

dlparkhurst:
Your rate is about 200 mol/s. Use PRINT statements to debug your RATES definition.

--- Code: ---#DATABASE C:\iCP\ORH Example - Silica Dissolution\phreeqc\LLNL_mod.dat
PHASES
MultiGlass
SiO2Al0.1326504481Fe0.0022Fe0.0257Mn0.0037772087Ca0.0027055058Na0.1864276569K0.0765685019(PO4)0.0001(OH)0.7316129321 + 0.7316129321H+ = SiO2(aq) + 0.1326504481Al+3 + 0.0022Fe+3 + 0.0257Fe+2 + 0.0037772087Mn+2 + 0.0027055058Ca+2 + 0.1864276569Na+ + 0.0765685019K+ + 0.0001PO4-3 + 0.7316129321H2O
log_k 6.6731
-delta_H -2.297 kJ/mol # Calculated enthalpy of reaction #from where
# Enthalpy of formation:
-analytic -623.7513 -0.0768 41744.3699 219.6997 -2695700.6140
# -Range: 0-300
RATES
MultiGlass
# parms(1): Aa = 4.03e-4 mol/m2/s
# parms(2): Ac = 0.105 mol/m2/s
# parms(3): Eaa = 45600 J/mol
# parms(4): Eac = 80000 J/mol
# parms(5): na = 0.30
# parms(6): nc = -0.4
# parms(7): R = 8.314 J/mol/K
# parms(8): sigma = 1
# parms(9): spec As = m2/kg
# parms(10): mass of MultiGlass = kg - for 1 mol = 85 g
-start
1 Aa = parm(1)
2 Ac = parm(2)
3 Eaa = parm(3)
4 Eac = parm(4)
4 na = parm(5)
6 nc = parm(6)
7 Rg = parm(7)
8 sigma = parm(8)

9 T = TK # convert ?C to K
10 Hydrogen = ACT("H+")

# pass specific surface in m?/kg:
11 sa = parm(9)

# Compute the saturation ratio for MultiGlass:
13 qok = SR("MultiGlass")

# Calculate the individual rate contributions:
14 k_a = Aa * (Hydrogen^na) * EXP(-Eaa/(Rg*T))
15 k_c = Ac * (Hydrogen^nc) * EXP(-Eac/(Rg*T))
16 rateplus = k_a + k_c

17 mass_rock = parm(10)
18 rate = rateplus * (1 - qok^(1/sigma)) * sa * mass_rock * 0.70 * GFW("SiO2") * 1e-3
#19 save rate
20 save rate * TIME
30 PRINT rate, sa, mass_rock, (1 - qok^(1/sigma))
-end
END
SOLUTION 0 #boundary water
temp 100
pH 11.0
units mol/kgw
density 1
Si 1e-6
Cl 1e-6
Na 0.1 charge
-water 1 # kg

SOLUTION 1-40
temp 100
pH 7
pe 4
redox pe
units mol/kgw
density 1
Na 1e-6
Cl 1e-6
Si 1e-6
-water 1 # kg

#EQUILIBRIUM_PHASES 1-40
# MultiGlass 0.0 3.5 dissolve_only

KINETICS 1-40
MultiGlass
#-formula
#-m
-m0 1
-parms 4.03e-4 0.105 45600 80000 0.30 -0.4 8.314 1 1e3 0.085
-time_step 1
-cvode
END

--- End code ---

Navigation

[0] Message Index

Go to full version