Model the water vapor fraction of a humid air inside a tower

Question

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.

clc

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);
    end

Ys



% Construct the desired plot
subplot(2,2,1)
plot(T1*Tr,Ys1);
xlabel('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Antonie Function');

subplot(2,2,2)
plot(T1*Tr,Ys2);
label('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Swietoslawski and Smith');

subplot(2,2,3)
plot(T1*Tr,Ys3);
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);
end

% 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);
end



% 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);
end