Implement an equilibrium solver

[jiweiqi]

This is very essential for a lot of the cases

[TJP-Karpowski]

Hi,
I am really interested in your project and would like to contribute to it:
so to get the equilibrium in Cantera most of the time we constrain "HP" and there the equilibrium is the minimum of Gibbs energy if I remember correctly.
Maybe one could use https://julianlsolvers.github.io/Optim.jl/stable/# for the minimisation of the Gibbs energy with the constrain of kepping the element fraction and the Enthalpy/pressure constant. What do you think? If I find the time I will look again into the CEA Online theory guide to refresh my memory on equilibrium calculations.

[TJP-Karpowski]

This is the current state of the equilibrium function. So while the result is not entirely wrong there is still this ln(X) term in the Gibbs energy which I currently neglect since it is nonlinear. I guess I have to look more closely at the paper you suggested to understand how they resolved that. I also tried to use the Nonlinear optimisation, but with little luck so far.

cd(@__DIR__)
using Arrhenius
using JuMP
using Clp
using GLPK
using Ipopt
using LinearAlgebra 
"""
determine the moles of each species in phase at equilibrium
given an array of initial moles X and a pressure and temperature
"""
function equilibriumTP(gas::Q,X0::V,temperature::B, pressure::B) where {Q<:Arrhenius.Solution,V<:AbstractArray,B<:Real}
    optimisation_problem = Model(Ipopt.Optimizer)
    @variable(optimisation_problem, 0<= X[1:gas.n_species] <= 1)
    N = sum(X0)
    X0 ./ N 
    enthaphy =get_H(gas, temperature, [], X0)
    entropy =get_S(gas, temperature, pressure, X0)
    elements0 = gas.ele_matrix * X0
    @constraint(optimisation_problem, gas.ele_matrix * X .== elements0)
    @constraint(optimisation_problem,  sum(X) == 1)
    
    #Definition of objective
    @expression(optimisation_problem, H, sum(X[i]*enthaphy[i] for i in 1:gas.n_species))
    @expression(optimisation_problem, S, sum(X[i]*entropy[i] for i in 1:gas.n_species))
    #@NLexpression(optimisation_problem, logX,  sum(X[i]*log(X[i]) for i in 1:gas.n_species))
    @expression(optimisation_problem, Gs, H - temperature * S )
    @objective(optimisation_problem, Min,  Gs ) # +logX
    display(optimisation_problem)
    optimize!(optimisation_problem)
    display(solution_summary(optimisation_problem))
    display(termination_status(optimisation_problem))
    return JuMP.value.(X)
end

gas = CreateSolution("../../mechanism/gri30.yaml")
ns = gas.n_species

T0 = 1200.0
P = one_atm 

X0 = zeros(ns);
X0[species_index(gas, "CH4")] = 0.5  
X0[species_index(gas, "O2")] = 1.0
X0[species_index(gas, "N2")] = 3.76
X0 = X0 ./ sum(X0);
Y0 = X2Y(gas, X0, dot(X0, gas.MW));

X_eq = equilibriumTP(gas,X0,T0,1e5) 
print(X_eq)```

[jiweiqi]

Your code looks nice to me. I am currently busy implementing a one-dimensional flame, more specifically gradient computing for one-D flame using auto-differentiation.

As for the reference, when I am looking for the one-d flame equations, I found that Chemkin uses the STAN lib, which is well documented at https://shepherd.caltech.edu/EDL/PublicResources/sdt/refs/STANJAN_1986_Reynolds.pdf.

[RSuryaNarayan]

Hello there, @TJP-Karpowski. This is Surya from India. I would like to know the status of the development of this feature and would be really interested in contributing. I was wondering if I should start from scratch following a re-implementation along the lines of @jiweiqi 's STANJAN reference or take you preliminary implementation using Optim.jl to a more complete one.

[TJP-Karpowski]

Hi @RSuryaNarayan,
so I have started with the Optim.jl and was a bit optimistic. It turns out that the problem optimizing the mole fractions directly is nonlinear due to the log(X) term in the entropy. Thats the reason for the Lagrangian Multiplier Method being used in STANJAN if I understood correctly. So I would say the best bet would be to re implement it based on the STANJAN reference. I also found this reference from NASAs CEA which describes a similar method: https://ntrs.nasa.gov/api/citations/19780009781/downloads/19780009781.pdf . So far I had no time to work on this feature further, so I appreciate your help :).

[RSuryaNarayan]

@TJP-Karpowski indeed. STANJAN uses an "element-potential" method to determine the mole-fractions based on the lagrange multipliers. I am yet to understand the document fully as I just started out. Regarding the other method using Optim.jl, I feel there would be some issues not only with the non-linear terms, but also with the initial guess. I've had little experience of coding an equilibrium solver for dissociating air (using simplified Lighthill's model) using NLsolve however the solution and convergence are quite sensitive to the initial guess. I am not sure whether this shall be mitigated completely using the method of element potentials, however, there is a mention of providing suitable initial guesses in the reference :)