Behavior of Peng Robinson equation of state for H-P queries

[Chris N]

I'm having an issue with making enthalpy-pressure queries using the Peng Robinson equation of state.

The setup of the problem is to fix the pressure at 101325 Pascals (the critical pressure is about 7MPa for reference) and select a low temperature of 100 Kelvin. At this state, use the setState_TP() method to set the state, which should be liquid. From that state query the enthalpy to get the starting enthalpy. Then use a setState_TP() query using a temperature of 500 Kelvin (critical temperature is 560 Kelvin). Repeat the previous process to get the ending enthalpy.  Then create an array of enthalpies to loop over and call the setState_HP() method using the enthalpy array entries.  From this I print out H,P,rho,T.


Within the ThermoPhase.cpp file there is the following method: setState_HPorUV
In the setState_HPorUV(double Htarget, double p, double rtol, bool doUV) method, for the HP option a strange thing appears to happen. During the Newton iteration there is a check on the convergence of the enthalpy, but that check is an OR conditional. And the other half of the conditional uses the temperature and the change in temperature as the metric for convergence. During cases where you're in the two phase region, the dt is essentially zero, and so that criteria erroneously signals a converged solution and the code continues on with a non-converged enthalpy. If I remove the temperature condition and keep the enthalpy condition, then the code properly tells me that the enthalpy has failed to converge.

Code from ThermoPhase.cpp:

// Convergence in H

double Herr = Htarget - Hnew;

double acpd = std::max(fabs(cpd), 1.0E-5);

double denom = std::max(fabs(Htarget), acpd * Tnew);

double HConvErr = fabs((Herr)/denom);

if (HConvErr < rtol || fabs(dt/Tnew) < rtol) { Is this a known issue? Am I doing something conceptually incorrect by making such queries? I've attached a sample code that runs the problem that I described.

[Ray Speth]

Hi Chris,

As we previously discussed, I don't think the current Peng-Robinson implementation actually represents two-phase mixtures. On top of that, though, I think you're right that the algorithm in ThermoPhase::setState_HPorUV wouldn't work for such a case anyway. This state setter is doing a nested iteration over temperature and pressure, and there's no way of setting of setting a two phase mixture to a particular value of the enthalpy when using those as your independent variables.

Regards,

Ray

[Chris N]

Thanks Ray. I'm looking through the MixtureFugacityTP class, and I notice that is some coding that uses the gibbs energy of the equation of state to determine the saturation pressure from a given temperature. It seems like maybe this class was going in that direction as some point?

The original algorithm in the ThermoPhase::setState_HPorUV takes in the P and H and first sets the density of the state using setPressure. Then it tries to find a temperature that generates the same enthalpy as the input value, H. And then sets this to be the temperature. This seems to work fairly fine. But it does struggle when even getting anywhere close to the two-phase dome(from the example the density starts of oscillate between vapor and gas even before getting to the saturation point).

I was thinking about making the following tweak by having the MixtureFugacityTP class override the setState_HPorUV() to use the following approach:
1.) Take Pressure from the input argument.
2.) Compute the saturation temperature at the given pressure.

3.) Using the pressure and saturation temperature, query the cubic eos to get the liquid and vapor enthalpies.

4.) Compare the input enthalpy to see if it lies between the liquid an vapor values, if so the system is two-phase.

5.) If the system is two-phase, then set the temperature to the saturation temperature.

6.) Use the liquid/vapor enthalpies to determine quality and then use that quality to estimate the density.

7.) Store the density of the mixture.
8.) If the enthalpy is outside of the bounds, then go back to the standard temperature iteration algorithm that was originally implemented.

Does this seem like a fundamentally unsound approach?

[Ray Speth]

Hi Chris,

I think in terms of the broad strokes this seems reasonable. One caveat is that you would need to implement setState_HP and setState_UV separately. From an implementation perspective, this is because only these methods are actually virtual methods that can be overridden in derived classes, and setState_HPorUV is a private method of the ThermoPhase class. But also, I think the procedure you outlined makes sense for the HP case but the UV case is going to be different.

I’d also suggest looking at the algorithms used by the Substance class (see src/tpx/Sub.cpp) that’s used for all the pure fluid equations of state, as I know those are fairly robust in the two-phase region.

If you’re able to get this working, I think this would make a valuable addition to these models in Cantera.

Regards,
Ray

​