Reactive Systems

[pw0908]

Started doing some work on trying to model reactive systems. Right now it will only work in a bulk system.

As an example:

components = [("ethanol", ["CH3" => 1, "CH2" => 1, "OH(P)" => 1]),
                          ("acetic acid", ["COOH" => 1, "CH3" => 1]),
                          ("ethyl acetate", ["CH3" => 1, "CH3COO" => 1, "CH2" => 1]),
                          ("water", ["H2O" => 1])]
reactions = [("ethanol"=>-1,"acetic acid"=>-1,"sodium"=>1,"chloride"=>1)]

eosmodel = UNIFAC(components)
model = ReactiveModel(["ethanol","acetic acid","ethyl acetate","water"], reactions; model = eosmodel)

equilibrate(model,1e5,298.15,[0.5,0.5,0.0,0.0]; z0 = [0.1])
4-element Vector{Float64}:
 0.4626300794049871
 0.4626300794049871
 0.03736992059501288
 0.03736992059501288

[pw0908]

Things that are left to do:

• Handle electrolyte reactive systems (i.e. aqueous systems). The reference system is different so a different activity needs to be used.
• Handle equilibrium systems. I'd imagine a 'cheat' would be to introduce equilibrate in-between steps?
• Create a database of formation energies.

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 155 lines in your changes missing coverage. Please review.

Project coverage is 77.36%. Comparing base (cf96c68) to head (e7b4161).
Report is 7 commits behind head on master.

Files	Patch %	Lines
src/models/Reactive/reactive.jl	0.00%	39 Missing ⚠️
src/models/Reactive/reactive_aq.jl	0.00%	39 Missing ⚠️
...c/methods/property_solvers/reactive/reactive_aq.jl	0.00%	29 Missing ⚠️
src/methods/property_solvers/reactive/reactive.jl	0.00%	23 Missing ⚠️
src/methods/pT.jl	0.00%	21 Missing ⚠️
src/models/Activity/equations.jl	0.00%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #284      +/-   ##
==========================================
- Coverage   77.93%   77.36%   -0.58%     
==========================================
  Files         316      321       +5     
  Lines       23585    23796     +211     
==========================================
+ Hits        18382    18410      +28     
- Misses       5203     5386     +183     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[pw0908]

Made some headway with the reactive aqueous systems. And got it to work to pure water to get pH:

species = ["water","hydronium","hydroxide"]
reactions = [("water"=>-2,"hydronium"=>1,"hydroxide"=>1)]

eosmodel = SAFTVREMie(["water"],["hydronium","hydroxide"])
model = ReactiveAqModel(species,reactions; model = eosmodel)

Clapeyron.pH(model,1e5, 298.15, [1.0, 0.0, 0.0]; z0=[1e-9]) # 7.0978252141976785

The difficult part here is the reference systems...

In the case of electrolytes, since the ions can't exist in isolation, the reference formation energies are measured at 1 mol/kg of the salt. This is fine when there is a single reaction (2 H2O<-> H3O+ + OH-) since that would just be the H2O 'salt'. However, if you have a weak acid, and thus two reactions (HA + H2O <-> H3O+ + A-), there are three ions, H3O+, A- and OH-. What's the reference system in this case? My feeling is I'll treat the reference as reaction-dependent. So the reference will change for each reaction, but at least the way the reference is treated will be consistent.

[longemen3000]

Maybe we could differentiate references by parameter name? (Gf_aq fvs Hf0) there are a finite amount of references anyway. I suppose we can convert between those if we have info of at least one of them

[pw0908]

So at the moment I'm differentiating them by the model type (ReactiveModel vs ReactiveAqModel). I think since the way the equations are solved and the concentration scales involved are different, it's best to keep it that way. However, I seem to have failed at defining the struct types since, ideally, ReactiveAqModel <: ReactiveModel, but I couldn't get it to work. Any advice @longemen3000?

[longemen3000]

i was just looking at the code, and the main difference between the aqueous and normal reactive model are the value of μ_ref. For example, activity_coefficient and activity can be defined as a function of μ_mix - μ_ref . i added an API that allows this explicit separation, maybe we could define transformations of references by modifying the value of μ_ref (at least for pure -> aqueous)