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

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 • 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. Dec 3, 2022 at 0:22
• 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. Dec 4, 2022 at 17:24