GetTWetBulbFromHumRatio does not expected results

[dmey]

I am using the GetTWetBulbFromHumRatio but noticed that the results are not as expected. when compared against EnergyPlus (PsyTwbFnTdbWPb, see this) or https://www.weather.gov/epz/wxcalc_rh, the results from GetTWetBulbFromHumRatio suggest an issue with the pressure.

E.g.

Temperature = 10.0 deg C
Humidity ratio = 0.0001 kg/kg
Pressure = [80000, 140000] Pa

 Pressure:    80000.000000000000
 GetTWetBulbFromHumRatio:   -1.5518645498801780      PsyTwbFnTdbWPb:  -0.88337304600021549
 --------------------------------------------------------------------
 Pressure:    90000.000000000000
 GetTWetBulbFromHumRatio:  -0.85689781113800689      PsyTwbFnTdbWPb:  -0.20562830125160234
 --------------------------------------------------------------------
 Pressure:    100000.00000000000
 GetTWetBulbFromHumRatio:   0.43708771369073951      PsyTwbFnTdbWPb:   0.40338928189695028
 --------------------------------------------------------------------
 Pressure:    110000.00000000000
 GetTWetBulbFromHumRatio:   0.98029150820405553      PsyTwbFnTdbWPb:   0.94753667228016969
 --------------------------------------------------------------------
 Pressure:    120000.00000000000
 GetTWetBulbFromHumRatio:    1.4604243562908890      PsyTwbFnTdbWPb:    1.4288966335071569
 --------------------------------------------------------------------
 Pressure:    130000.00000000000
 GetTWetBulbFromHumRatio:    1.8890212226606971      PsyTwbFnTdbWPb:    1.8584618711636722
 --------------------------------------------------------------------
 Pressure:    140000.00000000000
 GetTWetBulbFromHumRatio:    2.2742148244443414      PsyTwbFnTdbWPb:    2.2446875369332200
 --------------------------------------------------------------------

I have not looked at this yet in detail as to the issue is with the implementation or with and issue/limitation in the ASHRAE formulae.

[didierthevenard]

I don't believe this to be a bug, but rather result from the use of different formulae in the respective codes.

In PsychroLib:
In function GetSatHumRatio:
SatVaporPres = GetSatVapPres(TDryBulb)
SatHumRatio = 0.621945 * SatVaporPres / (Pressure - SatVaporPres)

In function GetHumRatioFromTWetBulb:
if TWetBulb >= 0:
HumRatio = ((2501. - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb))
/ (2501. + 1.86 * TDryBulb - 4.186 * TWetBulb)
else:
HumRatio = ((2830. - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb))
/ (2830. + 1.86 * TDryBulb - 2.1 * TWetBulb)

In EnergyPlus:
In function PsyTwbFnTdbWPb_raw:

Wstar=0.62198d0*PSatstar/(Patm-PSatstar)

! Calculate new humidity ratio and determine difference from known
! humidity ratio which is wStar calculated earlier
newW=((2501.0d0-2.381d0WBT)Wstar-(TDB-WBT)) / &
(2501.0d0+1.805d0TDB-4.186d0WBT)

The coefficients in PsychroLib are those used in the ASHRAE Handbook - Fundamentals over the past few editions. The coefficients in EnergyPlus are those found in Engineering Thermodynamics: SI Edition
By D.C. Look, G. Alexander, see here: https://books.google.ca/books?id=PI3pCAAAQBAJ&pg=PA448&lpg=PA448&dq=2501+2.381+1.805&source=bl&ots=ySfdCnLWh5&sig=ACfU3U2lelROZYvAsvafiG8cPgfe1o8_pw&hl=en&sa=X&ved=2ahUKEwjd_PKCirrgAhUixVkKHe8YAswQ6AEwAHoECAMQAQ#v=onepage&q=2501%202.381%201.805&f=false

I don't know if one set of coefficients is more accurate than the other, in any case PsychroLib is implementing ASHRAE and in that respect is probably if not correct, at least consistent.

[didierthevenard]

Further comments on this issue. The equations used by ASHRAE have changed over the years, for example the equation for enthalpy in the 2001 Handbook (which you can find here: http://www.ce.utexas.edu/prof/Novoselac/classes/ARE383/Handouts/F01_06SI.pdf) is:
h = 1.006t + W(2501+1.805t)
In the 2005 Handbook it became:
h = 1.006t+W(2501+1.86t)
Not sure where the jump in the coefficient of t came from.
By the looks of it, EnergyPlus and Look and Alexander are using the 2001 ASHRAE formulation, and furthermore assuming that the enthalpy of dry air is h_da = t whereas ASHRAE uses h_da= 1.006t.

[didierthevenard]

Last comment - the formulae for enthalpy mentioned above are actually an approximation of the data published in Table 3 of the psychrometrics chapter. A comparison of the 2001 and 2005 formulae to data from the table makes me wonder why the change was made, as the new formula appears less accurate than the old one. See attached file.
SatWaterVaporEnthalpyCalculation.xlsx

[dmey]

Perhaps this is not an issue per se but an inconsistency with other functions due to different coefficients or constants used, assumptions made, or validity checks. I will break this down into two points:

1. Coefficients and constants
   I think that the specific heat of water vapour is more in line to what I have commonly used in the past and seen. 1.86 kJ / kg/ K seem/ed reasonable to me but I have looked at the spreadsheet and noted that even though the error is smaller at values below 0 using the 'new' coefficient, the error appears to be larger at common climatic conditions. (e.g. 10-40). I have not investigated this in any detail but I would like to have a look at this in more details if I find the time in the next few days -- this is just a quick response.

Continuing from #16, I note that a number of constants that we use can be easily calculate from a set of more fundamentals constants (I will add a table of those shortly). The only issue with this, would be that the values would be slightly different than what is currently reported in the ASHRAE handbook. The question would therefore be whether we want 100% consistency or better accuracy and perhaps fewer empirical solutions when analytical solutions are available... What's your take on this?

2. I have noticed that under certain conditions our functions did not behave as expected therefore I tightened the validity checks. I have added these to a new branch for now. See this for the diff. If you are happy I will propagate the change to the other languages and merge this in -- this has improved some of the results when comparing against EnergyPlus already.

As a final comment for now, changing the value of specific heat of water vapour in PsychroLib will not significantly improved the results of the GetTWetBulbFromHumRatio function when tested under a large range of different conditioning. This may also be because in EnergyPlus, the function has added features/ different checked which we did not consider yet -- I note for example the calculation of boiling point of water (see this). I am not saying that we should implement those as well but see if there is an issue with any of our functions.

[dmey]

I have had a look at this issue in more details, the issue appears be purely numerical and, specifically to using the secant method rather than the false position method as done in EnergyPlus. There appears cases where we do not get convergence as we would do using the false position method (see this for example) -- this is probably why they chose to false position in EnergyPlus for this function.

I have created some tests for this function and attached the results (see psychrolib-v-eplus.zip). The findings can be summarised by the following histogram.

All significant errors are in the first few degrees -- this seems to indicates a convergence problem as given by using the bisection method. I have already tried in Fortran and we could simply change the loop and add a false position function that we can call. What do you think? This should be really straightforward if we do similarly to EnergyPlus but let me know if you had something else in mind we could use instead.
Re the other points I made in mu previous comments, I think we can continue with defining the constants directly in #16 to keep this issue tidy whilst I would like to propagate the other small changes about tightening function checks and to the other languages and merge it -- given though that these ceases are the same maybe I can add a helper function to take care of that to avoid repeating too much code... What do you think?

[didierthevenard]

Thanks for the various comments:

1. Tightening the validity checks: the code as proposed:
   do while ((abs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE) .or. (index < MAX_ITER_COUNT))
   will force the code to do MAX_ITER_COUNT loops at a minimum every time it does the calculation, even if convergence has been reached. Is that the idea?
2. Non-convergence. I am a bit confused where you are looking - whether it's in GetTWetBulbFromHumRatio (which uses a bisection method) or somewhere else? Can you provide a specific example of conditions for which convergence is not achieved? When you say 'under certain conditions our functions did not behave as expected' can you be more specific and provide an example?
3. Constants: I am okay with defining constants as long as we keep compatibility with ASHRAE formulae. The stumbling block and the one that in my mind you should focus on is whether constants can easily be defined in IP (I know they can be in SI).
4. All other changes: this library is an implementation of the ASHRAE formulae, It is designed to match what is in the ASHRAE Handbook, not what is in EnergyPlus, so I am not worried if the results from PsychroLib don't match what is in EnergyPlus. That said for higher resolution we should develop another library based for example on the IAPWS formulations.

[dmey]

Some quick comments:

1. There's a typo there -- it should be MIN_ITER_COUNT but yes, that should be the min number of loops to do in the iteration even if convergence is reached. That was not the only change in the branch (see https://github.com/psychrometrics/psychrolib/compare/dmey/tighten-validity-checks to compare the changes I propose) it also includes changes to bounding the humidity ratio value. I would suggest making this into a helper function so that I can simply be called when needed to avoid repetitions -- if agree, I am happy to implement this in all languages and create the PR.

2. That's correct -- I am suggesting to change the following (e.g. by using the false position method done in EnergyPlus) or something else which would give us convergence under any conditions.
   

   psychrolib/src/fortran/psychrolib.f90

   Lines 530 to 551 in 6f42ff1

   	TWetBulbSup = TDryBulb
   	TWetBulbInf = TDewPoint
   	TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.0
   	index = 1 ! initialise index
   	! Bisection loop
   	do while((TWetBulbSup - TWetBulbInf > PSYCHROLIB_TOLERANCE) .or. (index < MAX_ITER_COUNT))
   	! Compute humidity ratio at temperature Tstar
   	Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure)
   	! Get new bounds
   	if (Wstar > BoundedHumRatio) then
   	TWetBulbSup = TWetBulb
   	else
   	TWetBulbInf = TWetBulb
   	end if
   	! New guess of wet bulb temperature
   	TWetBulb = (TWetBulbSup + TWetBulbInf) / 2.0
   	index = index + 1
   	end do

   
   An example could be something like Temp = 7, Pressure = 130000, RH = 0.03
   You can also download the following test outputs I prepared which will summarise all the data and you can sort by max difference and easily find the conditions where convergence is not reached by the current implementation (root).

3. Good point -- I am not too familiar with Imperial units but let me have a look if I can find something.

4. Sure, it's just in case we find a discrepancy with a formula in the ASHRAE Handbook or find an obvious issue with it (e.g. maybe GetTWetBulbFromHumRatio) -- we already do that with GetTDewPointFromVapPres where we use ch. 1 eqn. 5 and 6 instead of eqn. 37 and 38 as they are much less accurate and have a narrower range of validity. Good point re the IAPWS formulations, though I don't think I'll have time to start working on those for another few months -- probably till the summer to be realistic.

[didierthevenard]

1. Before we do this I want to understand why it is necessary, so hold on while I convince myself of that.
   2 I have trouble understanding the 'something like Temp = 7, Pressure = 130000, RH = 0.03. Temp is the temperature, presumably in C. Pressure is the atmospheric pressure, presumably in Pa, but is that on Venus? It seems high. As for RH, is that relative humidity (RH) or humidity ratio (HR)?
2. Yes, IAPWS is down the road but not now. I am already devoting more time o this than I want to so I would like to keep changes etc to a minimum - only fix bugs.

[dmey]

1. OK
2. In the file I porvided you'll find all the values tested across a very large range (both inputs and outputs) -- not just the ones on Venus. But yes, Temperature [deg C], Pressure [Pa], Relative Humidity [1].
3. Sure

[didierthevenard]

Great thanks - but I still don't see the issue.

1. In the file you provided, it appears to me that you compare PsychroLib to EnergyPlus values. The changes to the code that you are suggesting are to solve a numerical issue - the iterative method currently used in PsychroLib not converging to the right solution, which is something different. What I want to see is an example of this: values of DBT, HR and Pressure that make GetTWetBulbFromHumRatio converge to the wrong solution.
2. You comparison file shows excellent agreement across all conditions for positive values of the wet bulb temperature (the max disagreement is 0.06 C which is completely acceptable). Where the larger disagreements happen is for negative values of temperature, where the disagreement can reach 0.7 C. But then again this is PsychroLib vs. EnergyPlus, I can't tell you which one is the most accurate. Does EnergyPlus have a case for negative temperatures? BTW in your spreadsheet you shouldn't calculate the error in relative terms since we are dealing with temperatures; absolute error is what should be used.

[dmey]

I think you hit the nail on the head! The error metric I used should have should have also included the absolute error or at lest used the relative error or relative percentage error scaled to Kelvin -- sorry! I agree with the peak at 0.7 degree C but increasing as we approach the positive values -- I have summarised this in the following scatter-plot. I have added a table below, were the errors are largest.

Temperature [°C]	Pressure [Pa]	Humidity Ratio [kg/kg]	TWB_psychrolib [K]	TWB_eplus [K]	Absolute error [K]	Relative error [1]	Relative Percentage Error [%]
12	80000	0.00001	272.4599218	273.1762074	0.716285582	0.002622064	0.26220643
12	80000	0.000109067	272.5963248	273.3097746	0.713449827	0.002610407	0.261040728
11	80000	0.00001	271.9696677	272.6664927	0.696825052	0.002555595	0.255559473
11	90000	9.07352E-05	272.7936297	273.4833545	0.68972483	0.002521999	0.25219993
11	80000	0.000102079	272.0989983	272.7855086	0.686510239	0.002516667	0.251666682
11	90000	0.000181497	272.9240433	273.611213	0.68716965	0.002511482	0.251148205
10	80000	0.00001	271.4699379	272.148061	0.678123105	0.002491743	0.249174329
11	80000	0.000204191	272.240863	272.9168316	0.675968592	0.00247683	0.247682999
10	80000	9.54882E-05	271.5916835	272.2606963	0.669012837	0.002457251	0.245725089
11	80000	0.000306337	272.3822682	273.0475064	0.665238195	0.002436346	0.243634598
9	80000	0.00001	270.9602157	271.6197503	0.659534582	0.002428154	0.2428154
10	90000	0.00001	272.1605965	272.8220732	0.661476707	0.002424572	0.242457181
10	80000	0.000191006	271.7273078	272.3859556	0.65864787	0.002418068	0.241806839
11	80000	0.000408517	272.5228425	273.178638	0.655795524	0.002400611	0.240061056
10	100000	7.63882E-05	272.8620563	273.5185834	0.656527036	0.002400301	0.240030139
10	90000	8.4877E-05	272.2707024	272.9238611	0.653158699	0.00239319	0.239319016
9	80000	8.92754E-05	271.0760204	271.7262255	0.650205075	0.002392868	0.239286831
11	80000	0.00051073	272.6627791	273.3162057	0.653426584	0.002390735	0.239073487
11	80000	0.000612977	272.8020524	273.4531904	0.651137974	0.002381168	0.238116795
10	80000	0.000286553	271.8623691	272.5106339	0.64826484	0.002378861	0.237886071
11	80000	0.000715257	272.9400469	273.5895953	0.649548396	0.002374171	0.237417068
8	80000	0.00001	270.4410331	271.0813598	0.640326686	0.00236212	0.236211994
9	90000	0.00001	271.6341351	272.2765882	0.642453089	0.002359561	0.235956052
9	80000	0.000178576	271.2047116	271.8456405	0.6409289	0.002357694	0.235769424
10	90000	0.000169777	272.3952503	273.0388002	0.6435499	0.002356991	0.235699065
10	80000	0.000382129	271.9964257	272.6347347	0.638308989	0.002341261	0.23412607
8	80000	8.34219E-05	270.5498163	271.181888	0.632071721	0.002330804	0.23308036
9	90000	7.93547E-05	271.7388498	272.3726781	0.633828364	0.002327063	0.232706294
9	80000	0.000267903	271.3334814	271.9645357	0.631054391	0.002320355	0.232035544
10	90000	0.0002547	272.5191308	273.1528735	0.633742698	0.002320103	0.232010262

Having said all this, I would agree in that one would not be able to tell you which one is the most accurate from this comparison... Is a disagreement of ~0.7 K significant for the application it is to be designed? In my specific case, I compared it against EnergyPlus as I started to notice large difference in my model outputs when switching from EnergyPlus to PsychroLib given that I needed a lightweight library to take care of the Psychometric calculations.

[didierthevenard]

I had a peek inside EnergyPlus but I don't really know how their ITERATE function works. In the body of the PsyTwbFnTdbWPb function though, it looks like they use only the formula appropriate for temperatures above zero, not an alternate formula for when temperatures are below zero as PsychroLib and ASHRAE do. Is that the source of the differences you observe between PsychroLib and EnergyPlus for negative temperatures?
In any case I am still not convinced we need to change anything to our code - I don't have any evidence that it is converging to the wrong value.

[dmey]

That's right, they only use the formula appropriate for temperatures above zero -- I tested the range again but this time using the formula for temperatures above zero for both below and above zero in PsychroLib and we get really good agreement against EnergyPlus. Max diff ~0.06 -- see below for plot. I agree, this is not a bug/issue with PsychroLib but actually a simplification in EnergyPlus. I have removed the bug label. Feel free to close this at your convenience.

[dmey]

For reference, I am including the tests output below.
psychrolib-hum-ratio-formula-above-zero-v-eplus.zip

[didierthevenard]

Issue is closed, no action needed.