MATLAB: For loop problem; could you help

forloop_variablesMATLAB

I wish to calculate w to plot for (T,w) for every thi values and T values, my goal is to plot psychrometric chart

do I need to use double for loop on this case??

eg for thi = 0,1 , T = 0:1:60 , find w ;; for thi = 0,2 , T = 0:1:60 find w2 ….. ;;; for thi = 1 , t = T = 0:1:60 find w(n)

pg(i) has been successfully computed but struggling on finding w vectors.

Could you please take a look into this?

%%% welxer formulation
g0 = -2.9912729*1000;
g1 = -6.0170128*1000;
g2 = 1.887643854*10;
g3 = -2.8354721*10^-2;
g4 = 1.7838301*10^-5;
g5 = -8.4150417*10^-10;
g6 = 4.4412543*10^-13;
g7 = 2.858487;
% T ; %%% 
%%%
T = 273.15:1:273.15+60 %%% temperature
thi = 0.1 : 0.1 : 1
   
    
for i =  1:length(T)
    
    Pg(i) = exp((g0/T(i)^2 + g1/T(i) + g2 + g3*T(i) + g4* T(i)^2 + g5 *T(i)^3 + g6*T(i)^4 ) + g7*log(T(i)))/1000;
                       %%%%% I would like save w vectors for for every thi = 0:0.1:1 , T = 0:1:60
                       w(i) = (0.622*thi*Pg(i)/101-thi*Pg(i);  %%%???
end

Best Answer

thi = 0.1:0.1:1;
T = 273.15:1:273.15+60; %%% temperature
N = numel(T);
w = nan(numel(thi),N); % preallocate!

Pg = nan(1,N);         % preallocate!
for k =  1:N
    Pg(k) = exp((g0/T(k)^2 + g1/T(k) + g2 + g3*T(k) + g4* T(k)^2 + g5 *T(k)^3 + g6*T(k)^4 ) + g7*log(T(k)))/1000;
    w(:,k) = 0.622*thi*Pg(k)./101-thi*Pg(k);
end

Every column of w is one result vector for one temperature.

https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html

Related Solutions

MATLAB: My fluid simulation exploding to very large numbers when I take small step values

I found a solution in the comments. I updated the value of n at each stage of the computations of j,k,l,o and now it is working as expected. Here is a code solution for the fluid simulation that works as intended:

clear;
clc;
format long g;
%% Defining the constants
L = 0.6; % Characteristic Length of evaporation surface in m
W = 0.6; % Width of the pan in m
A = L*W; % Surface area of one pan in m^2
d = 0.032; % Height of air columns in m 
rho = 1.225; % Density of Air in kg/m3
u = 5; % air speed in m/s
Ambient_Temp = 30;
Water_Temp = 40;
RH = 40;
Water_FlowRate = 6; % in LPH
Water_FlowRate = Water_FlowRate/3600; % Conversion to kg/s
%Calculating Absolute Humidity for w_g()
abs = (6.112 * exp((17.67 * Ambient_Temp)/(Ambient_Temp + 243.5)) * RH * 2.1674)/((Ambient_Temp+273.15)) * (1/1000);
ma = rho * u * W * d; % mass flow rate of the air (per pan) in kg/s
Kel = 273.15; % Celsius to Kelvin conversion
patm = 101325; % atmospheric pressure in Pa
delX = 0.01; % Size of each step
max_iters = 1/delX; % Number of Iterations
cpv = 1864; % cp of vapor in J/kg-K
cpa = 1005; % cp of dry air in J/kg-K
cpl = 4182; % cp of water in J/kg-K
ifgw = 2500 * 1000; % Latent heat of water at 0C in J/kg
% Constants in saturation pressure definition
a1 = -5800.2206;
a2 = 1.3914993;
a3 = -0.048640239;
a4 = 0.41764768 * 10^-4;
a5 = -0.14452093 * 10^-7;
a6 = 6.5459673;
%% Arrays
w_g = 1:max_iters;
T_g = 1:max_iters;
T_l = 1:max_iters;
ml = 1:max_iters;
ig = 1:max_iters;
psat = 1:max_iters;
iv = 1:max_iters;
w_int = 1:max_iters;
h = 1:max_iters;
Lef = 1:max_iters;
cpg = 1:max_iters;
k = 1:max_iters;
w_g(1) = abs/rho; % in relation to absolute humidity calculated above (w_g) is dimensionless so one must divide it by rho so the kg/kg cancels out
T_g(1) = Ambient_Temp + Kel; % Temperature of the Gas
T_l(1) = Water_Temp + Kel; % Temperature of the liquid
ml(1) = Water_FlowRate; % Mass flow rate of the liquid in Kg/s
ig(1) = (cpa)*(T_g(1) - Kel) + (w_g(1) * (ifgw + (cpv * (T_g(1) - Kel)))); % Enthalpy of the moist air in Joules/Kg
%% Dependent terms
va = ((0.009392 * T_l(1)) - 1.237) * 10^-5; % dynamic viscosity




%
Dh = (2 * W * d)/(W + d); % Hydraulic Diameter 
Re = (u * Dh)/ va; % Reynolds number calculation




Pr = (-0.0002697 * T_l(1)) + 0.81; %Prandtl number




n = 0.4;
lambda = 0.027; % Conductivity of air in W/m-K
psat(1) = exp((a1/T_l(1)) + (a2) + (a3 * T_l(1)) + (a4 * (T_l(1))^2) + (a5 * (T_l(1))^3) + (a6 * log(T_l(1)))); % Saturation pressure of the water

iv(1) = ifgw + cpv*(T_l(1) - Kel); % enthalpy of the water vapour

w_int(1) = (0.622 * psat(1))/(patm-psat(1)); % Moisture content at the interface

if (T_l(1) < T_g(1))
    n = 0.3;
end
h(1) = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh); % Defining the constant h

Lef(1) = ((0.865)^(2/3))*(((w_int(1) + 0.622)/(w_g(1) + 0.622) - 1)/(log((w_int(1) + 0.622)/(w_g(1) + 0.622)))); % Defining the Lewis constant

cpg(1) = cpa + w_g(1) * cpv; % Defining the specific heat of moist air 

k(1) = h(1)/(cpg(1) * Lef(1)); % Definition of constant K
%% Runge-Kutta 4th order Method 
for ctr = 2:max_iters
    
    
    % Defining first terms
    c = ctr - 1;
    
    ki1 = k(c);
    wint1 = w_int(c);
    wg1 = w_g(c);
    Tl1 = T_l(c);
    Tg1 = T_g(c);
    iv1 = iv(c);
    h1 = h(c);
    ml1 = ml(c);
    ig1 = ig(c);
    
    pd = wint1 - wg1;
    td = Tl1 - Tg1;
    
    j1 = delX * (-1 * ki1 * A * pd);
    k1 = delX * (ki1 * A * pd / ma);
    l1 = delX * ((h1 * A * td / ma) + (iv1 * ki1 * A * pd / ma));
    o1 = delX * ((A / (ml1 * cpl)) * ((ki1 * pd * (cpl*(Tl1 - Kel) - iv1)) - (h1 * td)));
    
    % Defining second terms
    
    ml2 = ml1 + (j1/2);
    wg2 = wg1 + (k1/2);
    ig2 = ig1 + (l1/2);
    Tl2 = Tl1 + (o1/2);
    
    va = ((0.009392 * Tl2) - 1.237) * 10^-5; % dynamic viscosity
    
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * Tl2) + 0.81; %Prandtl number
    
    psat2 = exp((a1/Tl2) + (a2) + (a3 * Tl2) + (a4 * (Tl2)^2) + (a5 * (Tl2)^3) + (a6 * log(Tl2)));
    iv2 = ifgw + cpv*(Tl2 - Kel);
    
    wint2 = (0.622 * psat2)/(patm-psat2);
    Tg2 = ((ig2 - (wg2*ifgw))/(cpa + wg2*cpv)) + Kel;
    
    if (Tl2 < Tg2)
        n = 0.3;
    end
    
    h2 = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh);
    
    Lef2 = ((0.865)^(2/3))*(((wint2 + 0.622)/(wg2 + 0.622) - 1)/(log((wint2 + 0.622)/(wg2 + 0.622))));
    cpg2 = cpa + wg2 * cpv;
    
    ki2 = h2/(cpg2 * Lef2);
    
    pd2 = wint2-wg2; 
    td2 = Tl2-Tg2;
    
    j2 = delX * (-1 * ki2 * A * pd2);
    k2 = delX * (ki2 * A * pd2 / ma);
    l2 = delX * ((h2 * A * td2 / ma) + (iv2 * ki2 * A * pd2 / ma));
    o2 = delX * ((A / ml2 / cpl) * (ki2 * pd2 * (cpl*(Tl2-Kel) - iv2) - (h2 * td2)));
    
    % Defining third terms
    
    ml3 = ml1 + (j2/2);
    wg3 = wg1 + (k2/2);
    ig3 = ig1 + (l2/2);
    Tl3 = Tl1 + (o2/2);
    
    va = ((0.009392 * Tl3) - 1.237) * 10^-5; % dynamic viscosity
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * Tl3) + 0.81; %Prandtl number
    
    psat3 = exp((a1/Tl3) + (a2) + (a3 * Tl3) + (a4 * (Tl3)^2) + (a5 * (Tl3)^3) + (a6 * log(Tl3)));
    iv3 = ifgw + cpv*(Tl3 - Kel);
    
    wint3 = (0.622 * psat3)/(patm-psat3);
    Tg3 = ((ig3 - (wg3*ifgw))/(cpa + wg3*cpv)) + Kel;
    
    if (Tl3 < Tg3)
        n = 0.3;
    end
    
    h3 = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh);
    
    Lef3 = ((0.865)^(2/3))*(((wint3 + 0.622)/(wg3 + 0.622) - 1)/(log((wint3 + 0.622)/(wg3 + 0.622))));
    cpg3 = cpa + wg3 * cpv;
    
    ki3 = h3/(cpg3 * Lef3);
    
    pd3 = wint3-wg3; 
    td3 = Tl3-Tg3;
    
    j3 = delX * (-1 * ki3 * A * pd3);
    k3 = delX * (ki3 * A * pd3 / ma);
    l3 = delX * ((h3 * A * td3 / ma) + (iv3 * ki3 * A * pd3 / ma));
    o3 = delX * ((A / ml3 / cpl) * (ki3 * pd3 * (cpl*(Tl3-Kel) - iv3) - (h3 * td3)));
    
    % Defining 4th terms
    
    ml4 = ml1 + (j3/2);
    wg4 = wg1 + (k3/2);
    ig4 = ig1 + (l3/2);
    Tl4 = Tl1 + (o3/2);
    
    va = ((0.009392 * Tl4) - 1.237) * 10^-5; % dynamic viscosity
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * Tl4) + 0.81; %Prandtl number
    
    psat4 = exp((a1/Tl4) + (a2) + (a3 * Tl4) + (a4 * (Tl4)^2) + (a5 * (Tl4)^3) + (a6 * log(Tl4)));
    iv4 = ifgw + cpv*(Tl4 - Kel);
    
    wint4 = (0.622 * psat4)/(patm-psat4);
    Tg4 = ((ig4 - (wg4*ifgw))/(cpa + wg4*cpv)) + Kel;
    
    if (Tl4 < Tg4)
        n = 0.3;
    end
    
    h4 = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh);
    
    Lef4 = ((0.865)^(2/3))*(((wint4 + 0.622)/(wg4 + 0.622) - 1)/(log((wint4 + 0.622)/(wg4 + 0.622))));
    cpg4 = cpa + wg4 * cpv;
    
    ki4 = h4/(cpg4 * Lef4);
    
    pd4 = wint4-wg4; 
    td4 = Tl4-Tg4;
    
    j4 = delX * (-1 * ki4 * A * pd4);
    k4 = delX * (ki4 * A * pd4 / ma);
    l4 = delX * ((h4 * A * td4 / ma) + (iv4 * ki4 * A * pd4 / ma));
    o4 = delX * ((A / ml4 / cpl) * (ki4 * pd4 * (cpl*(Tl4-Kel) - iv4) - (h4 * td)));
    
    % Finding the value at i+1 
    
    ml(ctr) = ml(c) + ((j1 + 2 * j2 + 2 * j3 + j4)/6);
    w_g(ctr) = w_g(c) + ((k1 + 2 * k2 + 2 * k3 + k4)/6);
    ig(ctr) = ig(c) + ((l1 + 2 * l2 + 2 * l3 + l4)/6);
    T_l(ctr) = T_l(c) + ((o1 + 2 * o2 + 2 * o3 + o4)/6);
    
    if (Tl1 <= 0) || (Tl2 <= 0) || (Tl3 <= 0) || (Tl4 <= 0)
        
        disp("ERROR!");
        disp("1 : " + wint1 + " psat : " + psat(c) + " Tl : " + Tl1 + " o value : " + o1 + " TL : " + Tl1);
        disp("2 : " + wint2 + " psat : " + psat2 + " Tl : " + Tl2 + " o value : " + o2 + " TL : " + Tl2);
        disp("3 : " + wint3 + " psat : " + psat3 + " Tl : " + Tl3 + " o value : " + o3 + " TL : " + Tl3);
        disp("4 : " + wint4 + " psat : " + psat4 + " Tl : " + Tl4 + " o value : " + o4 + " TL : " + Tl4);
        
        disp(ctr);
        break;
    end
    
   % Updating dependent values 
   
    T_g(ctr) = ((ig(ctr) - (w_g(ctr)*ifgw))/(cpa + w_g(ctr)*cpv)) + Kel;
    
    va = ((0.009392 * T_l(ctr)) - 1.237) * 10^-5; % dynamic viscosity
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * T_l(ctr)) + 0.81; %Prandtl number
    
    psat(ctr) = exp((a1/T_l(ctr)) + (a2) + (a3 * T_l(ctr)) + (a4 * (T_l(ctr))^2) + (a5 * (T_l(ctr))^3) + (a6 * log(T_l(ctr)))); % Saturation pressure of the water
    iv(ctr) = ifgw + cpv*(T_l(ctr) - Kel); % enthalpy of the water vapour
    w_int(ctr) = (0.622 * psat(ctr))/(patm-psat(ctr)); % Moisture content at the interface
    if (T_l(ctr) < T_g(ctr))
    n = 0.3;
    end
    h(ctr) = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh); % Defining the constant h
    Lef(ctr) = ((0.865)^(2/3))*(((w_int(ctr) + 0.622)/(w_g(ctr) + 0.622) - 1)/(log((w_int(ctr) + 0.622)/(w_g(ctr) + 0.622)))); % Defining the Lewis constant
    cpg(ctr) = cpa + w_g(ctr) * cpv; % Defining the specific heat of moist air 
    k(ctr) = h(ctr)/(cpg(ctr) * Lef(ctr)); % Definition of constant K   
   
end
%%
L_array = L/max_iters:L/max_iters:L;
t_w = 1/10000;
t = L/ma;
diff = ml(1) - ml(end);
evap = diff * 3600;
dw = w_int - w_g;
O_AH = w_g(end)* 1.225;
O_Tg = T_g(end) - Kel;
O_RH = (1000 * O_AH * (273.15 + O_Tg)) / (6.112 * exp((17.67 * O_Tg)/(O_Tg + 243.5)) * 2.1674);
disp("Output RH is : " + O_RH + "%");

MATLAB: Potting graphs with if else and for

If it throws an error, please include that information.

If that error is what I'm seeing:

Error using plot. Vectors must be the same length.

then the issue is that you're conditions don't include an x-value when t=4 (or -2,2, but these will be filled in with zeros). You can use <=, >=, or you don't want those points plotted, insert NaNs:

t = -4:0.1:4;
for k1 =1:length(t)
    if (t(k1)>-4) && (t(k1)<-2)
        x(k1) = (t(k1).^2) - (2.*t(k1))+3;
    elseif (t(k1)>-2) && (t(k1)<2)
        x(k1) = 4*cos(t(k1)*(2*pi*t(k1) - pi./8)) + 3*sin(2*pi*t(k1));
    elseif (t(k1)>2) && (t(k1)<4)
        x(k1) = sin(t(k1))./t(k1);
    else
        x(k1) = nan;
    end
end
plot(t,x)

Best Answer

Related Solutions

MATLAB: My fluid simulation exploding to very large numbers when I take small step values

MATLAB: Potting graphs with if else and for

Related Question