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

Question

I've been using this piece of code to solve an evaporation simulation, I tried using Euler's method to numerically solve the differential equations and this yielded results that were generally inaccurate (around 30% off from existing data). Still, with that method I got reasonably good results. To improve my calculations I tried using the Runge-Kutta method, but I kept getting answers in the 100,000s for my rate of evaporation (it should be less than 1 most of the time). Reducing the step value (stp) to around 0.1, I got pretty decent results; but they weren't that much more accurate than what I had got by using Euler's method (with Euler's method I used stp values of around 10^-6). When I used small values for stp, many of my terms turned complex, which didn't make sense.

I know this is a fairly long piece of code to sift through but I cannot figure out why this simulation explodes at small step values (e.g. below 0.01). I think it might be matlab rounding off the values at very small scales that might be throwing it off, but I can't be sure. Could anyone suggest methods that I could use to make this work. I think it might be an issue with matlab, but I wouldn't be surprised if there was a glaring issue with my code.

Thanks a lot in advance

Edit: I have added a small if statement below that I used as a flag (I added it after I noticed all my values turning complex). You can remove that if you want to run the entire simulation.

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 = 25;
Water_Temp = 40;
RH = 50;
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

stp = 0.001; % Size of each step

max_iters = 1/stp; % Number of Iterations

cpv = 1864; % cp of vapor in J/kg-K

cpa = 1000; % 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;
%




va = 1.6 * 10^-5; % dynamic viscosity

%
Dh = (2 * W * d)/(W + d); % Hydraulic Diameter 

Re = (u * Dh)/ va; % Reynolds number calculation

Pr = 0.705; %Prandtl number

n = 0.4;
lambda = 0.027; % Conductivity of air in W/m-K

%% 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

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
    
    c = ctr - 1;
    
    ki = k(c);
    wint = w_int(c);
    wg = w_g(c);
    Tl = T_l(c);
    Tg = T_g(c);
    iva = iv(c);
    hi = h(c);
    mli = ml(c);
    igi = ig(c);
    
    pd = wint - wg;
    td = Tl - Tg;
    
    % Calculation of j1 ,k1,l1,o1
    
    j1 = stp * (-1 * ki * A * pd);
    k1 = stp * (ki * A * pd / ma);
    l1 = stp * ((hi * A * td / ma) + (iva * ki * A * pd / ma));
    o1 = stp * (A / mli / cpl) * (ki * pd * (cpl*(Tl-Kel)-iva) - hi * td);
    
    % Calculation of j2,k2,l2,o2
    
    m2l = mli + (j1/2);
    w2g = wg + (k1/2);
    i2g = igi + (l1/2);
    T2l = Tl + (o1/2);
    
    
    p2sat = exp((a1/T2l) + (a2) + (a3 * T2l) + (a4 * (T2l)^2) + (a5 * (T2l)^3) + (a6 * log(T2l)));
    w2int = (0.622 * p2sat)/(patm-p2sat);
    
    Lef2 = ((0.865)^(2/3))*(((w2int + 0.622)/(w2g + 0.622) - 1)/(log((w2int + 0.622)/(w2g + 0.622))));
    cpg2 = cpa + w2g*cpv;
    ki2 = hi/(Lef2 * cpg2);
    
    T2g = ((i2g - (w2g*ifgw))/(cpa + w2g*cpv)) + Kel;
    i2v = ifgw + cpv*(T2l-Kel);
    
    p2d = w2int - w2g;
    t2d = T2l - T2g;
    
    j2 = stp * (-1 * ki2 * A * p2d);
    k2 = stp * (ki2 * A * p2d / ma);
    l2 = stp * (hi * A * t2d /ma) + (i2v * ki2 * A * p2d / ma);
    o2 = stp * (A / m2l / cpl) * (ki2 * p2d * (cpl * (T2l - Kel) - i2v) - (hi * t2d));
    
    % Calculation of j3,k3,l3,o3
    
    m3l = m2l + (j2/2);
    w3g = w2g + (k2/2);
    i3g = i2g + (l2/2);
    T3l = T2l + (o2/2);
    
    p3sat = exp((a1/T3l) + (a2) + (a3 * T3l) + (a4 * (T3l)^2) + (a5 * (T3l)^3) + (a6 * log(T3l)));
    w3int = (0.622 * p3sat)/(patm-p3sat);
    
    Lef3 = ((0.865)^(2/3))*(((w3int + 0.622)/(w3g + 0.622) - 1)/(log((w3int + 0.622)/(w3g + 0.622))));
    cpg3 = cpa + w3g*cpv;
    ki3 = hi/(Lef3 * cpg3);
    
    T3g = ((i3g - (w3g*ifgw))/(cpa + w3g*cpv)) + Kel;
    i3v = ifgw + cpv*(T3l-Kel);
    
    p3d = w3int - w3g;
    t3d = T3l - T3g;
    
    j3 = stp * (-1 * ki3 * A * p3d);
    k3 = stp * (ki3 * A * p3d / ma);
    l3 = stp * (hi * A * t3d /ma) + (i3v * ki3 * A * p3d / ma);
    o3 = stp * (A / m3l / cpl) * (ki3 * p3d * (cpl * (T3l - Kel) - i3v) - (hi * t3d));
    
    % Calculation of j4,k4,l4,o4
    
    m4l = m3l + (j3/2);
    w4g = w3g + (k3/2);
    i4g = i3g + (l3/2);
    T4l = T3l + (o3/2);
    
    p4sat = exp((a1/T4l) + (a2) + (a3 * T4l) + (a4 * (T4l)^2) + (a5 * (T4l)^3) + (a6 * log(T4l)));
    w4int = (0.622 * p4sat)/(patm-p4sat);
    
    Lef4 = ((0.865)^(2/3))*(((w4int + 0.622)/(w4g + 0.622) - 1)/(log((w4int + 0.622)/(w4g + 0.622))));
    cpg4 = cpa + w4g*cpv;
    ki4 = hi/(Lef4 * cpg4);
    
    T4g = ((i4g - (w4g*ifgw))/(cpa + w4g*cpv)) + Kel;
    i4v = ifgw + cpv*(T4l-Kel);
    
    p4d = w4int - w4g;
    t4d = T4l - T4g;
    
    j4 = stp * (-1 * ki4 * A * p4d);
    k4 = stp * (ki4 * A * p4d / ma);
    l4 = stp * (hi * A * t4d /ma) + (i4v * ki4 * A * p4d / ma);
    o4 = stp * (A / m4l / cpl) * (ki4 * p4d * (cpl * (T4l - Kel) - i4v) - (hi * t4d));
    
    % Plugging deriving next term 
    
    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;
    
    % Updating dependent terms

    
    T_g(ctr) = ((ig(ctr) - (w_g(ctr)*ifgw))/(cpa + w_g(ctr)*cpv)) + Kel;
     
    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
    
    if (imag(ml(ctr))) ~= 0 
        disp("ERROR!");
        break; 
    end
    
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;
plot(L_array, dw);

For reference, this is the code with Euler's method, it is much simpler, but continually underestimates the evaporation rate

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 = 25;
Water_Temp = 60;
RH = 50;
Brine_FlowRate = 6; % in LPH
Brine_FlowRate = Brine_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

stp = (10^-6); % Size of each step
max_iters = 1/stp; % Number of Iterations

cpv = 1864; % cp of vapor in J/kg-K
cpa = 1000; % 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;

%
va = 1.6 * 10^-5; % dynamic viscosity
%

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

lambda = 0.027; % Conductivity of air in W/m-K


%% 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) = Brine_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

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
%% Solving the DEs using Euler's method 

for ctr = 2:max_iters
    c = ctr - 1;
    
    ki = k(c);
    wint = w_int(c);
    wg = w_g(c);
    Tl = T_l(c);
    Tg = T_g(c);
    iva = iv(c);
    hi = h(c);
    mli = ml(c);
    
    pd = wint - wg;
    td = Tl - Tg;
    
    delM = stp * (-1 * ki * A * pd);
    delW = stp * (ki * A * pd / ma);
    delI = stp * ((hi * A * td / ma) + (iva * ki * A * pd / ma));
    delT = stp * (A / mli / cpl) * (ki * pd * (cpl*(Tl-Kel)-iva) - hi * td);
    
    ml(ctr) = ml(c) + delM;
            

            
            
            

    

        

	
    

        
Best Answer
    

    

								

								
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 + "%");

						    
						
		
    
        

    

Related Solutions
MATLAB: Energy balance issue with Simscape Heat Exchanger model

								

								
UPDATE: I resolved the issue with the energy balance. The heat transfer rate should be calculated as:
  Q  == epsilon * min(C1,C2) * (T_A1 - T_A2);
And also, I was using counter flow heat exchanger effectiveness formula instead of parallel flow heat exchanger effectiveness. That caused the unphyscial temperature profile of the fluids at outlet. The system and code are both working as expected now. I validated the code and system with the library component of Heat Exchanger (TL-TL). 
Hope this helps anyone trying to do similar thing!
						    
						
MATLAB: For loop problem; could you help

								

								
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 Question

• How can get the value for each time step for a variable inside of a function
• Basic for loop problem
• My function will not work to calculate the ‘Cumene inlet temperature’ section

Answer 1

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 + "%");