Estimation of psychrometric properties using Mathematica ​
​Authors: Housam Binous, André Von-Held Soares and Ahmed Bellagi

We compute psychrometric properties using Eqs. (1)–(11) of Ref. 1 and compare our data versus their counterparts published in Refs. 1-2 (see Table 2 reproduced below) for five case studies.​​
​​In order to estimate psychrometric properties, only three input parameters are required. These are usually the dry bulb temperature, dbt, and total pressure, pt, as well as one of the following parameters: (a) wet-bulb temperature, wbt, (b) dew-point temperature, tdp, (c) relative humidity, rh, (d) partial pressure of water-vapor, pw, (e) absolute humidity, w, (f) humid specific heat, cpma, and (g) specific volume, va. ​Although a few authors use a sequential approach for the determination of the unknown psychrometric properties. Such approach can be somewhat tricky for the novice. Here, we adopt a much simpler and general strategy by solving a set of nonlinear algebraic equation (Eqs. 1-11 of Ref. 1) using the built-in Mathematica command : FindRoot. Hence, all unknown properties are determined simultaneously in just one unique calculation independently of the adopted three input parameters. ​References :​[1] Gupta, P . K . and Patel, R . N . (2017), A teaching–learning tool for elementary psychrometric processes on psychrometric chart using MATLAB . Comput Appl Eng Educ, 25 : 458 - 467. https://doi.org/10.1002/cae.21813[2] Beltrán-Prieto, J.C., Beltrán-Prieto, L.A. and Nguyen, L.H.B.S. (2016), Estimation of psychrometric parameters of vapor water mixtures in air. Comput Appl Eng Educ, 24: 39-43. https://doi.org/10.1002/cae.21670​
In[]:=
eq1[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=cpma==1.005+1.88w;​​​​eq2[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=va==0.287(273.16+dbt)/(pt-pw);​​​​eq3[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=h==(1.006+1.775w)dbt+2501w;​​​​eq4[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=mu==100w/ws;​​​​eq5[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=w==0.622pw/(pt-pw);​​​​eq6[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=pw==rhpws/100;​​​​eq7[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=ws==0.622pws/(pt-pws);​​​​eq8[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=wbt==dbtArcTan[0.151977(rh+8.313659)^0.5]+ArcTan[dbt+rh]-4.686035-ArcTan[rh-1.676331]+0.00391838rh^1.5ArcTan[0.023101rh];​​​​eq9[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=tdp==Evaluate[If[-50<=tdp<=0,5.994+12.41Log[pw]+0.4273Log[pw]^2,If[0<=tdp<=50,6.983+14.38Log[pw]+1.079Log[pw]^2,​​13.80+9.478Log[pw]+1.991Log[pw]^2]]];​​​​eq10[dbt_,pt_,pw_,pws_,w_,rh_,tdp_,wbt_,cpma_,va_,ws_,h_,mu_]:=pws==Evaluate[If[-40<=dbt<=0,Exp[24.2779-6238.64/(dbt+273.16)-0.344438Log[dbt+273.16]],Exp[-7511.52/(dbt+273.16)​​+89.63121+0.023998970(dbt+273.16)​​-1.165451110^-5(dbt+273.16)^2​​-1.281033610^-8(dbt+273.16)^3​​+2.099840510^-11(dbt+273.16)^4​​-12.150799Log[dbt+273.16]]]];
CASE 1: given are the dry bulb temperature, the total pressure and the partial pressure of water-vapor
In[]:=
myfun1[dbt_,pt_,pw_]:=Module[{sys,unknowns},sys={eq1[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq2[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq3[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq4[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq5[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq6[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq7[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq8[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq9[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq10[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu]};​​unknowns={{pws,3},{w,0.03},{rh,50},{tdp,20},{wbt,20},{cpma,0.75},{ws,0.03},{h,100},{mu,0.9},{va,0.85}};​​FindRoot[sys,unknowns]]
In[]:=
CASE1=myfun1[26.7,101.325,2.76]
Out[]=
{pws3.50453,w0.0174171,rh78.7553,tdp22.6941,wbt23.7844,cpma1.03774,ws0.0222838,h71.2459,mu78.1604,va0.873128}
CASE 2 given are the dry bulb temperature, the total pressure and absolute humidity
In[]:=
myfun2[dbt_,pt_,w_]:=Module[{sys,unknowns},sys={eq1[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq2[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq3[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq4[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq5[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq6[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq7[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq8[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq9[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq10[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu]};​​unknowns={{pws,3},{pw,3},{rh,50},{tdp,20},{wbt,20},{cpma,0.75},{ws,0.03},{h,100},{mu,0.9},{va,0.85}};​​FindRoot[sys,unknowns]]
In[]:=
CASE2=myfun2[87.8,101.325,0.03]
Out[]=
{pws64.4776,pw4.66219,rh7.23072,tdp31.6781,wbt42.8907,cpma1.0614,ws1.08841,h168.032,mu2.75632,va1.07172}
CASE 3 given are the dry bulb temperature, the total pressure and relative humidity
In[]:=
myfun3[dbt_,pt_,rh_]:=Module[{sys,unknowns},sys={eq1[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq2[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq3[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq4[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq5[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq6[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq7[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq8[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq9[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq10[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu]};​​unknowns={{pws,4},{w,0.013},{pw,2},{tdp,18},{wbt,22},{cpma,0.75},{ws,0.022},{h,70},{mu,0.9},{va,0.85}};​​FindRoot[sys,unknowns]]
In[]:=
CASE3=myfun3[30.0,101.325,50.0]
FindRoot
:The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.
Out[]=
{pws4.24532,w0.0133091,pw2.12266,tdp18.4177,wbt22.2968,cpma1.03002,ws0.0272002,h64.1748,mu48.9301,va0.877065}
CASE 4 given are the dry bulb temperature, the total pressure and the dew-point temperature
In[]:=
myfun4[dbt_,pt_,tdp_]:=Module[{sys,unknowns},sys={eq1[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq2[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq3[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq4[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq5[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq6[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq7[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq8[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq9[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu],eq10[dbt,pt,pw,pws,w,rh,tdp,wbt,cpma,va,ws,h,mu]};​​unknowns={{pws,3},{w,0.03},{rh,50},{pw,2},{wbt,20},{cpma,0.75},{ws,0.03},{h,100},{mu,0.9},{va,0.85}};​​FindRoot[sys,unknowns]]
CASE 5 given are the dry bulb temperature, the total pressure and the wet-bulb temperature