Calculation not succeeded, but why?

[Damfow]

I am doing some calculations of solution of acid-alkaline gases in acid-base aqueous solution. A simple HCl(g) in air vs NaOH(aq) solution is simple and the calculation is defined as "Succeeded", but I want to go further.

The problem that I wanted to solve is H2S(g) + NH3(g) + CO2(g) vs NaOH(aq). I use the speciation method to know the involved components and my own list of components, but both solution are defined as "Not succeeded", but I do not know the reason, so I am handcuffed with no clue of how to do.

My scope is to figure out myself why is not solved. Is it possible to get more info about the conditions that the solver consider as "Not succeeded"?

[allanleal]

Hi @Damfow - the issue you are facing could be caused by many reasons (e.g., use of inadequate thermodynamic database and activity model for the problem, input conditions, etc.). Without a minimal reproducible example, it's difficult to assist you (but it seems you want to figure out by yourself).

The solver fails to compute the solution if, by default, the calculation does not converge in 200 iterations (i.e., the residuals of all linear and non-linear equations governing the equilibrium state of the system do not decrease by 1e-8 by default). You can check the log of the calculation using:

options = EquilibriumOptions()
options.optima.output.active = True

# Using EquilibriumSolver
solver = EquilibriumSolver(system)  # or using specs
solver.setOptions(options)  

result = solver.solve(state)

# Using equilibrate
result = equilibrate(state, options)

You'll get a new local text file named optima.log.txt by default.

[Damfow]

Hi again, thank you for your super fast answer.

Now, obviously, I need some legend to understand the data that pulls out.

I see some output parameters with names: error, ex_max, ep_max, ew_max, x[-], y[-], s[-], ex[-], ey[-]. I guess that those with initial "e" correspond to "error". So these should be less than 1e-8, right?

Anyway, I think it is better to ask you for a right thermodynamic package for my system... Are supcrtbl or phreeqc correct for my system?

[allanleal]

Yes, those values starting with e correspond to (residual) error. ex[H+] is the residual of the first-order optimality equation corresponding to H+ (i.e., the equilibrium equation for H+).

Is supcrtbl or phreeqc correct for my system?

It's hard to answer this because I don't know the thermodynamic conditions under which you're doing your calculations. Every database is different; they are usually optimized against a certain number of experimental data. For your gas phase, we may even need to define suitable binary interaction parameters for the gases (for increased accuracy).

Maybe your choice of database is not the cause of your problem; it could probably be how you are formulating the problem with Reaktoro. Sharing how you are doing this (with a simpler/simulated problem, if the original problem statement is confidential) will allow us to advise you better. This will also help us understand how the interface of Reaktoro could be made more intuitive.

[Damfow]

It is a simple (or not...) acid-base equilibrium that happens inside an scrubber (liquid-gas contactor). Temperature and pressure are almost ambient (slightly heated).

Here some piece of the code:

db = SupcrtDatabase("supcrtbl")
# db = PhreeqcDatabase("phreeqc.dat")
# db = PhreeqcDatabase("pitzer.dat")

# %% Initial values of the gas stream
Temp = 60 # ºC
V_gas = 12_500 # m3/h @T
VN_gas = V_gas * 273 / (273 + Temp) # Nm3/h
Conc_H2O = 29.8 # %vol
Conc_CO2 = 30.2 # %vol
Conc_O2 = 0.0001 # %vol
Conc_H2S = 0.5 # %vol
Conc_NH3 = 3.0 # %vol
Conc_N2 =  100-Conc_H2O-Conc_CO2-Conc_O2-Conc_H2S-Conc_NH3 # %vol

n_H2O = Conc_H2O / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_CO2 = Conc_CO2 / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_O2 = Conc_O2 / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_H2S = Conc_H2S / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_NH3 = Conc_NH3 / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_N2 = Conc_N2 / 100 * VN_gas / (273.15 * 0.082) # kmol/h

# Scrubber conditions
Q_H2O = 440 # kg/h water
Conc_NaOH = 20 # %w solution 
m_H2O = Q_H2O * (1-Conc_NaOH/100) # kg/h
m_NaOH = Q_H2O * Conc_NaOH / 100 # kg/h

Temp_Sc = 30 # ºC, Temp. scrubber

aq_comps = "H2O(aq) H+ OH- Na+ NaOH(aq) HCO3- CO3-2 CO2(aq) "
aq_comps += "H2S(aq) HS- "
aq_comps += "NH3(aq) NH4+"

gas_comps = "CO2(g) H2O(g) N2(g) H2S(g) NH3(g)"

solution = AqueousPhase(aq_comps)
gas = GaseousPhase(gas_comps)
system = ChemicalSystem(db, solution, gas)

state = ChemicalState(system)

solver = EquilibriumSolver(specs)
aprops = AqueousProps(state)

state.setTemperature(Temp_Sc, "celsius")
state.setPressure(1.0, "atm")

state.set("H2O(aq)", m_H2O, "kg")
state.set("NaOH(aq)", m_NaOH, "kg")
state.set("H2S(aq)", n_H2S, "kmol")
state.set("NH3(aq)", n_NH3, "kmol")

state.set("H2O(g)", n_H2O, "kmol")
state.set("CO2(g)", n_CO2, "kmol")
state.set("N2(g)", n_N2, "kmol")

# Solving...
# result = equilibrate(state)
result = solver.solve(state)
# result = solver.solve(state, conditions)

print("EQ STATE")
print(state)
aprops = AqueousProps(state)

print("AQUEOUS PROPERTIES AT EQUILIBRIUM")
print(aprops)

print("pH: ", aprops.pH())
print("Succeesful?", result.succeeded())

[allanleal]

Thanks for sharing the code above. From the log of the calculation I could identify the issue. It has to do with how the residual errors are checked and the scale of the species amounts you are using. The calculation essentially converged, but the largest residual could not be decreased below 1.88098e-08 (the default tolerance is 1.0e-8).

The code below is a revised version of yours to resolve the convergence check issue (tolerance increased to 1e-6).

I'll investigate more this and make the necessary changes in Reaktoro. I'll also check how Reaktoro v2.8.2 (latest version with some improvements in Optima) behaves. I used Reaktoro v2.7.2 for testing your code.

from reaktoro import *

db = SupcrtDatabase("supcrtbl")

# %% Initial values of the gas stream
Temp = 60 # ºC
V_gas = 12_500 # m3/h @T
VN_gas = V_gas * 273 / (273 + Temp) # Nm3/h
Conc_H2O = 29.8 # %vol
Conc_CO2 = 30.2 # %vol
Conc_O2 = 0.0001 # %vol
Conc_H2S = 0.5 # %vol
Conc_NH3 = 3.0 # %vol
Conc_N2 =  100-Conc_H2O-Conc_CO2-Conc_O2-Conc_H2S-Conc_NH3 # %vol

n_H2O = Conc_H2O / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_CO2 = Conc_CO2 / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_O2 = Conc_O2 / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_H2S = Conc_H2S / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_NH3 = Conc_NH3 / 100 * VN_gas / (273.15 * 0.082) # kmol/h
n_N2 = Conc_N2 / 100 * VN_gas / (273.15 * 0.082) # kmol/h

# Scrubber conditions
Q_H2O = 440 # kg/h water
Conc_NaOH = 20 # %w solution 
m_H2O = Q_H2O * (1-Conc_NaOH/100) # kg/h
m_NaOH = Q_H2O * Conc_NaOH / 100 # kg/h

Temp_Sc = 30 # ºC, Temp. scrubber

aq_comps = "H2O(aq) H+ OH- Na+ NaOH(aq) HCO3- CO3-2 CO2(aq) "
aq_comps += "H2S(aq) HS- "
aq_comps += "NH3(aq) NH4+"

gas_comps = "CO2(g) H2O(g) N2(g) H2S(g) NH3(g)"

solution = AqueousPhase(aq_comps)
gas = GaseousPhase(gas_comps)
system = ChemicalSystem(db, solution, gas)

state = ChemicalState(system)

solver = EquilibriumSolver(system)

state.setTemperature(Temp_Sc, "celsius")
state.setPressure(1.0, "atm")

state.set("H2O(aq)", m_H2O, "kg")
state.set("NaOH(aq)", m_NaOH, "kg")
state.set("H2S(aq)", n_H2S, "kmol")
state.set("NH3(aq)", n_NH3, "kmol")

state.set("H2O(g)", n_H2O, "kmol")
state.set("CO2(g)", n_CO2, "kmol")
state.set("N2(g)", n_N2, "kmol")

options = EquilibriumOptions()
options.optima.output.active = True
options.optima.convergence.tolerance = 1e-6
solver.setOptions(options)

result = solver.solve(state)

assert result.succeeded()

print("EQ STATE")
print(state)
aprops = AqueousProps(state)

print("AQUEOUS PROPERTIES AT EQUILIBRIUM")
print(aprops)

print("pH: ", aprops.pH())
print("Succeesful?", result.succeeded())

Please note that accuracy is a different story. Your setup above uses ideal activity models for the aqueous and gaseous phases (i.e., activity coefficients are 1 for aqueous species, and fugacity coefficients are 1 for the gases).

Let me know if you encounter further issues.

[allanleal]

In the meantime, I checked how Reaktoro v2.8.2 behaves for the problem above. It reduces the error to 3.8667e-09 and thus no need to tweak the tolerance value.

Note: After executing conda update reaktoro, I noticed that an existing dependency, spdlog, kept its version, and an error happended about libspdlog1.12 could not be found. I had to do conda install spdlog=1.12.

[Damfow]

I was working with a version of Reaktoro 2.6. Now my caclulations are precise. Thank you so much.

Anyway, I think that my problem could help other users to check if the final answer is near or far from the solution or it is totally worng.

Is there any reference to a example for checking the different variables of the parameters obtained from the variable result?

[allanleal]

Thanks for reporting that the newer version is working well for your problem above. Let me know if you find newer issues.

For your question above, I did not understand what you mean by "different variables of the parameters".

[Damfow]

The variable "result" have some parameters that are useful for checking the solution. I would like to have some example of how to access to every parameters of "result". It will be great if there is a webpage with all the parameters listed.

[allanleal]

Got it. result is of type EquilibriumResult. Here is the API Reference for this class:

https://reaktoro.org/api/structReaktoro_1_1EquilibriumResult.html#a39fce64e37d8a2493f4ff764708d958a

Right now, not much is exposed there from optima member, which has type Optima::Result and its docs need to be merged with that of Reaktoro.