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

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;
    ig(ctr) = ig(c) + delI;
    T_l(ctr) = T_l(c) + delT;
    w_g(ctr) = w_g(c) + delW;
    
    T_g(ctr) = ((ig(c) - (wg*ifgw))/(cpa + wg*cpv)) + Kel;
    
    % Updating dependent terms
    
    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;
plot(L_array, dw);

[Cris LaPierre]

When ctr = 323, p4sat is greater than patm. This results in w4int being negative. As part of your calculation of Lef4, you take the natural log of w4int (plus a very small value). The natural log of a negative number is complex.

[Anirudh Bhalekar]

Thanks a lot for the insight. I am looking into the code and trying to redo the calculations for the constants j,k,l,o. The weird part is that with high stp numbers (like.1 or even .05) it doesnt do anything too drastic, but with smaller step numbers, thhe terms become complex.

My initial guess was that matlab was internally rounding off calculations that resulted in values that were too small which resulted in compounding errors that made the arguments for the log function negative, but that assumption is contradicted by the fact that I can run the Euler's method code just fine with extremely low stp values (think 10^-5 and orders of magnitude lower).

[Anirudh Bhalekar]

Alright, I found the mistake. I ran the check for T_l < T_g in each of my computations of the constants in the RK method and accordingly updated h. It works just fine now. Thanks a lot for pointing out the negative log error, it really helped me pinpoint which section of the code was messing up.

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 = 40;
Water_Temp = 60;
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

delX = 0.01; % Size of each step
max_iters = 1/delX; % 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
    
    
    % 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);
    
    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);
    
    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);
    
    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 next value 
    
    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;
    
    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;
plot(L_array, dw);

Answer 1

[Anirudh Bhalekar]

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