Gibbs Reactor Calculation Methods

[Alex M]

I was going through the Gibbs Reactor tutorial, and I stumbled across an apparent issue.

I've modeled the reaction in the tutorial using the ChemSep components with both the Raoult's Law and UNIFAC thermo models. Both give the same answers. I've attached the Raoult's Law version since it runs so much faster. I've also modeled it with the DWSIM components (other than carbon monoxide, which is from CheResources in the tutorial, but I don't have that database so I used ChemSep again) and gotten very similar answers (off by a small percentage). (Is there an easy way to get the CheResources database without extensive repurposing of Pedro Fajardo's Excel Add-In?)

Regardless of the components used or the thermo model (provided the thermo model works), I get a discrepancy between the "DirectMinimization" and "ReactionExtents" calculation methods. The outlet streams come out essentially identical (even the same specific enthalpy), but the Heat Load is significantly different. I've got 410 kW with the DirectMinimization method (which is very close to the number in the tutorial), but 34 kW with the ReactionExtents method. Try as I might, I've been unable to determine a reason for the difference (they match in the tutorial).

It's worth noting that multiplying the mass flow by the difference in specific enthalpies (outlet-inlet) yields the 34 kW number the ReactionExtents method calculates (within 0.2 kW).

Have I done something wrong in my configuration?

Three other observations:

1. The Peng-Robinson EOS does not work as used in the version 1.7 tutorial (unless there is something really special about that CheResources carbon monoxide data). If you have the carbon monoxide present, it puts all the components in the vapor outlet stream but in the liquid phase. Even just creating a Material Stream with those conditions yields a stream in the liquid phase (not what I expected at 1000 K).

2. When I save and close the simulation then reopen it later, at some point the Gibbs Reactor in the ReactionExtents method becomes unconfigurable. I have to delete it and insert a new one. I've tried to save it, close and immediately attach this time so hopefully it won't be messed up yet.

3. The ReactionExtents method does not seem to work if you've not run the DirectMinimization method first. The error says the destination array is not long enough to copy all the items in the collection. Is this intentional?

[Daniel Medeiros]

Fixed the reaction extents method, both methods now result in 410 kW.

About the Carbon Dioxide thing, it seems that it is missing vapor pressure information in the ChemSep database, I'll take a look.

If the reactors become unconfigurable, try calculating them first, then edit.

There was a bug in the reaction extents method, the variable that stores the number of reactants wasn't being updated, it was always set to 0.

The next build will contain the fixes listed above.

Regards,
Daniel

[Alex M]

Can you tell me how you're calculating the Initial vs. Final Gibbs Free Energy (shown in the Gibbs Reactor block after the run is over), and how you're calculating the heat load?

I expected (using d for delta) dG = dH - TdS. Using the specific enthalpies and specific entropies tabulated for the inlet and outlet streams, I'm not coming up with the difference in Gibbs Energy indicated in the reactor block.

Along the same lines, I expected the heat load calculated by the isothermal reactor to equal the difference in enthalpies (Hreaction = Hproducts - Hreactants), but again using the specific enthalpies of the inlet and outlet streams I get a different value than that indicated by the Gibbs reactor block (and it's affiliated energy stream).

I do see that the direct minimization and reaction extents methods are getting the same answer now, which is good. I can also run the reaction extents method without doing the direct minimization method first, so well done there.

The PR EOS still thinks the outlet stream is liquid.

[Daniel Medeiros]

Gibbs energy is being calculated as the sum of the ideal gas gibbs energy of formation (DGf = DHf - T*DSf) and the non-ideal part through the fugacity coefficients (mole sum).

Heat load = Hp - Hr + DHr

Hp = products enthalpy
Hr = reactants enthalpy
DHr = heat of reaction

DHr is calculated from the sum of the ideal gas enthalpy of formation multiplied by the mole number variation of each component. Hp and Hr already include the ideal gas part. If you choose the reaction extents method, it will be calculated from the heat of reaction multiplied by the reaction extent value.

If you choose the adiabatic option, the heat of reaction will be added to the products enthalpy and the temperature will change as a result.

You can double check the procedures directly in-code, look for the Calculate() method: https://github.com/DanWBR/dwsim3/blob/master/DWSIM/Objects/Reactors/Gibbs.vb

Maybe the phase identification algorithm is identifying the phase as being liquid-like, take a look here: http://dwsim.inforside.com.br/wiki/index.php?title=Phase_Identification_Parameter

[Alex M]

Alright, I've managed to generally satisfy myself on the heat load calculation.

And you're right about the Phase Identification Parameter. If I turn that setting off the PR EOS produces a vapor, essentially matching the NRTL, UNIFAC and Raoult's Law results.

Is there any reason a reaction would run with the Reaction Extents method but not Direct Minimization? I've been playing with a methanol synthesis (CO+2H2=CH3OH), and I've observed that putting it in direct minimization produces an error, which says it "reached a stationary point of the objective function (singular gradient matrix)" and doesn't solve.

I attached the simulation.

On another note, I did manage to almost exactly match literature data for a catalyzed reaction at those conditions using the equilibrium reactor and the literature's Keq value. Obviously I wouldn't expect a Gibbs reactor to achieve that when there's a catalyst affecting the conversion.

[Daniel Medeiros]

Hi Alex,

There isn't a direct reason for that except for the convergence method which is completely different. The direct minimization technique uses lagrange multipliers and works only with a vapor phase, maybe that's why you're not getting anything on your methanol synthesis.

I did some debugging of your simulation and saw that the initialization using ideal chemical potentials doesn`t show any methanol formation at the specified conditions - hence the error. Maybe this is related to the absence of a catalyst. In this case you'd have to heat the mixture to compensate, but that may not work also.

I remember someone asking for catalyst support, I already have the literature references but hadn`t the time to do that yet, unfortunately, but it is on the TODO list.

Regards,
Daniel