Adiabatic Flame temperature computation

[Sush1090]

Hello @longemen3000 and @pw0908

I would like to know that based on all the current implementations in Clapeyron.jl can there be an extension towards computing adiabatic flame temperature? By any chance does there already exist an implementation somewhere? I am not sure if Catalyst.jl already has one.

My question arises from the implementation in Cantera.
Here is their working example: link.

My understanding of the underlying thermodynamics is not of the highest order hence I post the question here first before implementing anything.
Do correct me if I am wrong

So can we do the following:

1. To do this we can provide two EoSModel's input (reactants) and output (products) with known mole fractions as we know the reaction.

2. Then we solve for

$$p(V_{products} - V_{reactants}) = U_{products} - U_{reactants}$$

This will probably be a two variable NL solve problem based over Clapeyron.VT_internal_energy and we solve for volume and temperature. But there is a small catch, that the enthalpy of products and reactants has to be same. Hence

$$ h_{products} = h_{reactants} $$

https://en.wikipedia.org/wiki/Adiabatic_flame_temperature

Let me know what you think? Is this sensible or am I wrong regarding the thermodynamics?

[longemen3000]

for constant volume adiabatic flame temp, that seems ok. The main problem is providing the compositions of reactants and products (because we don't have a reaction framework), but if the user provides it, then we can just add a solver for that

[Sush1090]

My idea is to provide only one EoSModel that contains both products and reactants. And pass z_prod and z_reactants as parameters. I am assuming the user knows the reaction.

[longemen3000]

fair enough, one possible problem is a procedure for finding the initial points, probably supposing ideal gas could help?

[Sush1090]

Yes I think we need to use ideal gas.

There is an implementation here :

https://github.com/kyleniemeyer/computational-thermo/blob/main/book%2Fcontent%2Fchemical-equilibrium%2Fadiabatic-flame-temperature.ipynb

I think this is what we need. We probably also need to pass elemental composition.

[Sush1090]

Hello @longemen3000 ,

I was trying to reproduce the implementation for this to Clapeyron from here:

https://github.com/kyleniemeyer/computational-thermo/blob/main/content/chemical-equilibrium/adiabatic-flame-temperature.ipynb

But I have not been able to find a convergent solution. There are generally two type of problems I end up in:

1. Generic one where solution does not converge
2. The iterator goes in domains internally where we end up taking sqrt of negative values.

I am not sure if my EoSModel for this application is correct.

Here is my code (I am trying to reproduce what was in the Jupiter notebook)

using Clapeyron

# H4 + 0.5O2 ->  H2O 
using Unitful
import Unitful: K, Pa, J, mol

const atm1 = 101325.0Pa

components = ["hydrogen","oxygen","Water"]
model = cPR(components,idealmodel = BasicIdeal)#,reference_state = :nbp)

pressure = 101325.0*10Pa # Pa
T_init = 1000.0K  # K

moles_init = [1.0, 2.0, 0.0].*mol  # initial moles of H2, O2, H2O

elemental_comp = [
     2 0 2 ;  # H
     0 2 1    # O
]

function lagrange_system(x,pressure,model::EoSModel,elemental_comp::Matrix,T_init,moles_init)
    @assert size(elemental_comp,1) == 2 " There are only products and reactants hence elemental composition matrix should have 2 rows"
    moles = [x[1],x[2],x[3]].*oneunit(1.0mol)
    multipliers = [x[4],x[5]]*oneunit(1.0J/mol)  # Lagrange multipliers for H and O balance
    temperature = x[6]*oneunit(1.0K)
    mole_fractions = moles ./ sum(moles)
    # pures = split_model(model)
    gibbs = zeros(length(model.components))*oneunit(1.0J/mol)
    h_init = zeros(length(model.components))*oneunit(1.0J/mol)
    h_final = zeros(length(model.components))*oneunit(1.0J/mol)
    chemical_potentials = zeros(length(model.components))
    
    for (idx,comp) in enumerate(components)
        z_ = zeros(length(model.components)).*oneunit(1.0mol) # change z internally to not mess with reference states 
        z_[idx] = 1.0mol
        gibbs[idx] = Clapeyron.gibbs_free_energy(model, atm1, temperature, z_,phase = :g)/oneunit(1.0mol) # gibbs per mole of component
        h_final[idx] = Clapeyron.enthalpy(model, atm1, temperature, z_,phase = :g)/oneunit(1.0mol)  # enthalpy per mole of component
        h_init[idx] = Clapeyron.enthalpy(model, atm1, T_init, z_,phase = :g)/oneunit(1.0mol)  # enthalpy per mole of component
    end
    # R = Clapeyron.Rgas()*(1.0J/(mol*K))
    chemical_potentials = Clapeyron.chemical_potential(model, pressure, temperature, moles,phase = :g)
    initial_mole_elements = elemental_comp * moles_init
    mole_elements = elemental_comp * moles
    init_enthalpy = sum(h_init .* moles_init)
    final_enthalpy = sum(h_final .* moles)
    # chemical_potentials = gibbs .+ R * temperature .* log.(mole_fractions*pressure/atm1)  # avoid log(0)


    res_moles = mole_elements .- initial_mole_elements
    res_enthalpy = final_enthalpy - init_enthalpy
    res_chem_pot = zeros(length(chemical_potentials))*oneunit(1.0J/mol)#chemical_potential_res(model,pressure,temperature,moles,phase = :g)#zeros(length(chemical_potentials))*oneunit(1.0J/mol)
    for i in eachindex(chemical_potentials)
        res_chem_pot[i] = chemical_potentials[i] + sum(elemental_comp[:,i] .* multipliers)
    end
    return vcat(ustrip.(res_moles),ustrip.(res_chem_pot),ustrip.(res_enthalpy))
end
f(x) = lagrange_system(x,pressure,model,elemental_comp,T_init,moles_init)

using NLsolve
x0 = [1.0, 1.0, 1.0, 100.0, 100.0, 2200.0]
sol = nlsolve(f, x0,method = :trust_region, factor = 0.1,autoscale = false)

Any tips on this would be really helpful.

Thank you