I am going to model the moisture air inside an energy tower, therefore. I formulated the following equation to calculate the fraction of water vapor, Y in the moist air, where $\rho_w$ is the water vapor density,$\rho$ is the total density $T$ is the temperature, $R_a, R_w$ are the gas constant in dry air aid water vapor respectively. And $p_w$ and $p^s_w$ are the water vapor pressure and saturated vapor pressure \begin{align} Y &= \frac{\tilde{\rho}_w}{\tilde{\rho}} = \frac{\tilde{p}_w}{R_w\tilde{\rho}\tilde{T}} = \frac{\tilde{p}_w}{R_w\tilde{T}(\tilde{\rho}_w+\tilde{\rho})} = \frac{\tilde{p}_w}{R_w\big(\frac{\tilde{p}_w}{R_w}+\frac{\tilde{p}_a}{R_a}\big)}= \frac{\tilde{p}_w}{p_w+\big(\tilde{p}-\tilde{p}_w\big)\frac{R_a}{R_w}} \nonumber \\ &= \frac{R_a}{R_w}\frac{\tilde{p}_w}{\big(\tilde{p}-\tilde{p}_w\big(1-\frac{R_a}{R_w}\big)\big)} = \frac{R_a}{R_w}\frac{\phi\cdot\tilde{p}^s_w}{\big(\tilde{p}-\phi\cdot\tilde{p}^s_w\big(1-\frac{R_a}{R_w}\big)\big)} \end{align}

with the relative humidity, $\phi$ \begin{equation*} \phi = \frac{\tilde{p}_w}{\tilde{p}^s_w} \end{equation*} Obtaining saturation function with relative humidity being unity. i.e. $\phi=1$ and after the pressure split for the small Mach number correction $p=p_0=1$.(not necessary in this question). I obtain the saturation function, $Y^s$ as follows \begin{equation}\label{1} Y^s(T) = = \frac{R_a}{R_w}\frac{ p^s_w(T)}{\big(1-\big(1-\frac{R_a}{R_w}\big)\cdot p^s_w(T)\big)} \end{equation}

So I expect the $Y^s$ to be 1 when the air is saturated with vapor. So I tried to model this with several saturated vapor pressure functions $p^s_w$. as follows

  1. Antonie Function\begin{equation} p^s_w(T) = ptr\cdot exp\Big(17.28-\frac{4102.99}{T-35.719}\Big) \end{equation} where $ptr=611.657 Pa$ is the triple point of the water.

  2. Emporical Vapor pressure equations, Swietoslawski and Smith (1938)\begin{align} log P_w^s(T) &= A - \frac{B}{T}+\frac{CX}{T}(10^{DX^2}-1)-E(10^{-FY^{5/4}}) \end{align} where $X=T^2 - 2.937x10^5; A= 8.3074650; B=2.0051x10^3; C=1.3869x10^{-4}; D = 1.1965x10^{-11}; E=4.4x10^{-3}; F = 5.7148x10^{-3}; Y= 647.27 - T$

  3. Saturation vapor pressure of water and Ice, J. Huang (2018) \begin{align} P^s_w = \frac{exp(34.494-\frac{4294.99}{T+237.1})}{(T+105)^{1.57}} \text{ for $T>0$} \end{align}

Following is the Matlab code which I used to calculate this $Y^s$ using the above saturation functions. But my problem is I did not receive the expected saturation curve here usually got negative values or a constant. The tests are done with the temperature data with $\tilde{T}\in[273.16,333.16]$ and reference temp $T_r=T(0,t)=296.5K$ as the temperature at the bottom of the tower. What I expect here is that $Y^s$ grows with the temperature $T$ and after certain temperature, $Y^s$ should be constant(ideally 1). But I don't get the required results. I don't know actually my code is wrong or I'm using inappropriate function.


Rw = 461; Ra = 287.5;  % gas constant of water wapor and dry air
Tr = 296.5;

Pws = Ra/Rw;
Ptr = 611.66; % triple point pressure of water

% Testing Saturation function

% Temperature range

 T2 = sort(randi([273 400],1,30));% generate a random temperature distribution
 T1 = T2/Tr;

n = size(T1,2);
x = linspace(0,1,n); 

Ys1 = zeros(1,n);
Ys2 = zeros(1,n);
Ys3 = zeros(1,n);

% saturation pressure calculation
    for j=1:n
      Ys1(j,n) = antonie(T1(j),Tr,Pws);
      Ys2(j,n) = ss(T1(j),Tr,Pws,Ptr);

      Ys3(j,n) = ijh(T1(j),Tr,Ptr,Pws);


% Construct the desired plot
xlabel('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Antonie Function');

label('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Swietoslawski and Smith');

xlabel('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Ice, J. Huang');

% saturation pressure functions

function g = antonie(T,Tr,pws)
    T= T*Tr;
    ptr = 611.657;
    g = ptr*exp(17.2799-(4102.99/(T-35.719)));
    Ys = pws*g/(1-(1-pws)*g);

% According to the paper: Derivation and use of Antonie function for the
% temp. range 0-374(C°)
function Ys = ss(T,Tr,pws,ptr)
    T = T*Tr;
    A = 8.3074650; B = 2005.1; C= 1.3869*10^(-4); D= 1.1965*10^(-11);
    E = 4.4*10^(-3); F = 5.7148*10^(-3);
    Y = 374.11-T; X = T^2-2.937*10^5;
    g = ptr*exp(A-(B/T)+(C*X/T)*(10^(D*(X^2))-1)-E*(10^(-F*(Y^(5/4)))));
    Ys = pws*g/(1-(1-pws)*g);

% Refer : https://journals.ametsoc.org/view/journals/apme/57/6/jamc-d-17-0334.1.xml?tab_body=pdf
function ys = ijh(T,Tr,ptr,pws)
    T = T*Tr;
    t = T+237.1;
    g = ptr*exp(34.494-(4924.99/t))/((T+105)^1.57);
    ys = pws*g/(1+(1-pws)*g);

matlab plots

  • 1
    $\begingroup$ You're unlikely going to get good answers here, for two reasons: (i) You don't actually explain what you expect and what you see. You're only saying 'can someone explain to me', from which I assume that you expect a certain behavior but get something else, but you should be explicit about this because we're not all experts in the humidity of air. (ii) You dump a lengthy piece of code on us, but what does it actually do? Are we supposed to look over it and find a bug? If that's your goal, you need to explain what it does and why it doesn't produce what you expect it to. $\endgroup$ Commented Dec 3, 2022 at 0:22
  • $\begingroup$ Notice that your question doesn't actually contain a question mark. Phrasing your question so that it does often helps clarify what you need for yourself and others. $\endgroup$
    – Richard
    Commented Dec 4, 2022 at 17:24


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.