DanWBR / dwsim

DWSIM is a Steady-State and Dynamic Sequential Modular Chemical Process Simulator for Windows, Linux and macOS.
https://dwsim.org
GNU General Public License v3.0
297 stars 96 forks source link

Heat capacity / enthalpy of hypothetical component #52

Closed andr1976 closed 3 years ago

andr1976 commented 3 years ago

I suspect something might be off with the estimation of polynomial coefficients for the calculation of heat capacity/temperature dependant properties.

I have tried adding a hypothetical component with NBP, MW, Pc, Tc, Zc and omega. I get significantly lower mass specific heat capacity 0.27 KJ/kg K vs 2.32 kJ/kg K also the Cp/Cv values seems high 1.45 vs 1.04 from other simulator/thermpackage.

P = 1 atm T = 200 C

The hypo is C7 like fraction.

MW = 108.47 Tc = 302.45 C Pc = 2789.28 kPa Accentric= 0.3265 Vc = 0.4470 m3/kmole

No group contribution methods selected

hypo1.zip

DanWBR commented 3 years ago

How did you generate this compound?

andr1976 commented 3 years ago

I used the compound creator wizard, but maybe this unintended for pseudo component?

DanWBR commented 3 years ago

yep. without entering UNIFAC structure information or Joback groups DWSIM has no way of estimating ideal gas Cp. When you use the Petroleum fraction characterization utilities, it calculates Watson K, which can be used by a Lee-Kesler correlation to estimate Cp.

andr1976 commented 3 years ago

Hmm, I am still far off with respect to density (DWSIM 6.3 kg/m3 vs 2.87 ) (was close before) and still off on Cp (DWSIM 1,25 vs 2.33 kJ/kg K, very strange since the Pure component Cp for ideal gas is quite close to 2.3 in DWSIM), compressibility now also off (DWSIM 0.87 vs 0.97). petrol_hypo.zip

DanWBR commented 3 years ago

You have some options to choose from

Captura de Tela 2020-11-24 às 18 30 07
andr1976 commented 3 years ago

I know we might be comparing oranges with apples here, but the diffs are still quite significant. We use "phase behavior of pertroleum reservoir fluids" Pedersen, Christensen for calculating coefs for Cp (4).
Yeah, but Pc, Tc, omega, SG, MW is quite close, so I suspect the main diff is the method for Cp - and somehow compressibility is strange

andr1976 commented 3 years ago

Also the ideal gas Cp looks fine, but is significantly below in the stream tab - I wouldn't expect such a profound real gas effect on Cp image

DanWBR commented 3 years ago

Perhaps there is something wrong in here?

https://github.com/DanWBR/dwsim6/blob/fae9175fa51c710d5b4151ca728937fc26137a52/DWSIM.Thermodynamics/Property%20Packages/Models/FluidProperties.vb#L100-L118

andr1976 commented 3 years ago

I can check, but Cp_ig seems ok. Ok, I am sorry, I must be tired, I just tried setting Tc, Pc methods to Lee-Kessler and now we are getting there, density much closer and Cp also much closer. It is 20 to midnight here in Denmark, so off to bed, I will look more into this tomorrow

andr1976 commented 3 years ago

Ok, fresh eyes on the problem. It seems to give reasonable results now for Cp, Density and compressibility. There seems to be an error in the critical pressure when estimated by Lee-Keesler by a factor of 10

Another thing: Speed of Sound is enormous! 3436-5000 m/s depending on property methods.

DanWBR commented 3 years ago

https://github.com/DanWBR/dwsim6/blob/e3c195b394c17208a082aa687121249df07f5fd1/DWSIM.Thermodynamics/Property%20Packages/PropertyPackage.vb#L836-L851

DanWBR commented 3 years ago

Bulk modulus is the inverse of Isothermal Compressibility, which depends on the selected model.

https://github.com/DanWBR/dwsim6/blob/e3c195b394c17208a082aa687121249df07f5fd1/DWSIM.Thermodynamics/Property%20Packages/PropertyPackage.vb#L755-L808

andr1976 commented 3 years ago

My quick and dirty handcalc gives me approx. 200 m/s based on SQRT(Cp/Cv R T/M), back calculating with a vapor density of 2.86 kg/m3 gives me an isothermal compressibility 2 orders magnitude lower than stated in the stream tab

Maybe Dim K As Double = 1 / P0 - 1 / Z (Z1 - Z) / (P1 - P0) shall be Dim K As Double = - 1 / Z (Z1 - Z) / (P1 - P0)

hypo_case.zip

andr1976 commented 3 years ago

is Z really Z (compressibility), or is it actually the molar volume returned from AUX_Z? the K formula is on volume, but your symbols might indicate compressibility factor?

DanWBR commented 3 years ago

image

DanWBR commented 3 years ago

https://petrowiki.org/Isothermal_compressibility_of_gases

andr1976 commented 3 years ago

Ahh ok, I just checked the basic definition. Also stated: For gases at low pressures, the second term is small, and the isothermal compressibility can be approximated by cg ≈ 1/p. i.e. 1/1.013e5 approx. 9e-6 1/Pa, DWSIM gives 1.37858E-08 1/Pa

DanWBR commented 3 years ago

Because it always calculates the second term?

You can get a report of the internals by enabling the solution inspector.

andr1976 commented 3 years ago

It was just to get a reality check of the estimated magnitude of the isothermal compressibility, which still seems to be off by 1.5-2 orders of magnitude compared to hand calcs. Real gas effects should be fairly minimal at the chosen conditions, so the handcalcs should be a good benchmark for a sanity check.

DanWBR commented 3 years ago

Try the solution inspector. It will be easier to spot the error since you have the correct values at hand.

andr1976 commented 3 years ago

Ok I will try

andr1976 commented 3 years ago

It doesn't get me any closer to finding the root cause - maybe next step will be running DWSIM from source :-) The isothermal compressibility is like a solid and not a gas

andr1976 commented 3 years ago

Everything looks right in the code - but my wildest guess would be that the wrong root of Z is returned actually corresponding to a liquid phase

DanWBR commented 3 years ago

Possibly. And you're right, the inspector report has nothing for isothermal compressibility 😅

DanWBR commented 3 years ago

Isothermal Compressibility is calculated analytically for PR/SRK. The error is probably in here:

https://github.com/DanWBR/dwsim6/blob/10f7f21a6b9390e4ef030686498e5af14c552bdc/DWSIM.Thermodynamics/Base%20Classes/MichelsenBase.vb#L43-L80

DanWBR commented 3 years ago

If we're not able to find the error I'll switch to the general numerical method.

andr1976 commented 3 years ago

Methane at RT looks sane, 392 m/s, 1/K = 9.9e-6 1/Pa

DanWBR commented 3 years ago

calculated numerically:

Captura de Tela 2020-11-25 às 10 21 55
andr1976 commented 3 years ago

That seems more like it - strange though that the analytically derived quantity isn't right. Just tried Octane at I get 4328 m/s at 200 C and ambient. Maybe the number for methane was on the low side also (i find a tabulated value of 446 m/s at RT) - so the suspicion might be that something is off in the analytical expression - i might look more into this another day.

andr1976 commented 3 years ago

Now for the J-T coefficient.... ahh just kidding, that will be another day :-)

DanWBR commented 3 years ago

I'll use numerical derivatives for now. It doesn't seem to be that much slower.

andr1976 commented 3 years ago

I wonder why there's no difference in the analytical expressions for PR and SRK? I found this in a thesis for the Joule-Thomson, but 1/(dV/dP) is also defined image image

andr1976 commented 3 years ago

I am closing this for now, since we now get expected values. The above can be used for inspiration, when revoking the analytical method for isothermal compressibility/J-T coefficient. I have made a separate case for the Lee-Kessler Method for Pc

andr1976 commented 3 years ago

I have tried to look a little more into the isothermal compressibility. I tried implemeting the above and also reinstanted the analytical methods in the MichelsenBase.vb both seems to give around the same as the numerical method. Hmm, which I think is a little strange given that I previously got much too high values with the analytical method.