How do I estimate average of a variable over the discretized domain

Dear Community, I have written th code (as attached) to compute water adsorption in silica gel bed. In the following part of the code I am calculating the relative humidity, phi.

for i = 1:n
    P_sat = 610.78*exp((17.27*(Tair(i)-273.15))/(Tair(i)-35.85)); 
    P_vap=(W(i)*R*rho_air*Tair(i))/M_w; 
    phi= (P_vap*100)/P_sat; % relative humidity
    Q_eq = ((a*b1*phi)/1+(b1*phi)) + ((a*b2*(phi^q))/(1+(b2*(phi^n)))); % Langmuir-sips isotherm equation 
    DQDt(i) = K_LDF*(Q_eq - Q(i));   %LDF
end

In the equation of P_vap; W is the absolute humidity and Tair is the air temperature. The current codes uses W(i) and Tair(i) at each specifc n, where n is the number of discretization. I have used n=100, i.e., I have divided the whole air channel (length 0.2mm) into 100 divisions. The attached excel sheets are the outout of W and Tair.

Instead of W(i) and Tair(i), I would like to use the average of W and Tair for the discretized domain (i.e., the average value of columns B, C, D, E...in the excel sheets) at each time step in the equation of P_vap. Similarly, in the equation of P_sat, I want to use the average of Tair for the domain discretized in 'n' divisions.

I also want to outout two other parameters:1. average of Q for the 'n' discretized domain integrated at each time step (i.e., mass of water adsorbed in the bed), and 2. difference of the average of W of the 'n' discretized domain at the inlet and outlet (i.e., dehumidification capacity).

How can I implement this in the above section of the code?

I shall be highly grateful if someone can help.

[Torsten]

Why do you want to do this ? It's wrong because the adsorption equilibrium is a local value - thus depends on the state in cell i, not on the average over all cells.

[JUBAIR]

I want to plot the total mass of water adsorbed (Q_total =mean of Q(i=1:n)) in the bed vs. time. How can I do that?

Also, how do I estimate the difference of W between inlet and outlet at each time step?

Thanks for your help.

[Torsten]

How do you get the mass of adsobed water in cell i from Q ? I guess the unit of Q is not kg, but kg/kg or mol/kg, isn't it ? If it's in kg/kg, then I would guess that the total mass of water adsorbed is sum_i (V_i*(1-epsilon)*rho_s*Q_i) where V_i is the volume of cell i, rho_s is density of adsorbens and epsilon is porosity of the bed.

[JUBAIR]

Hi Toresten, you are right. Unit of Q is kg/kg. Let me try to use your formula. I will get back to you if I am successful. Thanks.

[Torsten]

Be careful with rho_s. Usually, one uses rho_s = m_s/(V_s + V_p) as the solid density in the simulation where V_s is the volume of the (pure) solid and V_p is the intra-particle volume. That's what I used above in computing the mass of adsorbed water.

Since you don't use the finite-volume method, I guess, the formula from above can be adjusted for your case as

(1-epsilon) * rho_s * ((V_1 * Q_1) / 2 + ( sum_{i=1}^{i=N-1} (V_i*Q_i + V_(i+1)*Q_(i+1)) / 2 ) + V_N*Q_N / 2)

Answer 1

[Marie Anna NOVIELLO]

• Calculate W and Tair averages: Rather of utilizing W(i) and Tair(i) for each division, compute the average value of W and Tair throughout the discretised domain at each time step.
• Use the following averages in the equations: Replace W(i) and Tair(i) with the derived averages when computing P_vap and P_sat to represent the overall domain conditions at each time step.
• Calculate the total mass of water adsorbed. Integrate (average) the Q values over the domain at each time step to calculate the total mass of water adsorbed in the bed.
• Calculate the dehumidification capacity: Calculate the difference between the average W at the inlet and outlet (first and last discretization points), which represents the dehumidification capacity.

% Initialize output arrays
Q_total = zeros(1, n); % to store total adsorbed water in the bed
W_dehumidification_capacity = zeros(1, n); % to store dehumidification capacity

for i = 1:n
    % Calculate average W and Tair for the discretized domain
    avg_W = mean(W_data(:, i)); % Average of W for the spatial domain at time step i
    avg_Tair = mean(Tair_data(:, i)); % Average of Tair for the spatial domain at time step i

    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85)); 
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w; 
    phi = (P_vap * 100) / P_sat; % relative humidity

    % Langmuir-sips isotherm equation
    Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    
    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));   

    % Calculate total mass of water adsorbed in the bed (integrated over the domain)
    Q_total(i) = sum(Q) / n; 

    % Calculate dehumidification capacity
    W_inlet = W_data(1, i); % Inlet (first column in W_data)
    W_outlet = W_data(n, i); % Outlet (last column in W_data)
    W_dehumidification_capacity(i) = mean(W_outlet) - mean(W_inlet); 
end

[Torsten]

@Marie Anna NOVIELLO

Do you mean "i" to loop over the time steps or the grid points ? If it's a loop over the output times, I don't think it's reasonable to use averaged values Q(i) and Q_eq to compute the adsorbed mass. After the simulation has finished, Q as a function of the spatial coordinate is given for each output time. So why not use these values directly to get the adsorbed mass ?

[JUBAIR]

@Marie Anna NOVIELLO I think your understanding is correct. We must use W and Tair average to compute RH. I have implemented your code in the main MATLAB file but getting some errors. I am atatching the m file and snapshot of error message herewith. In addition, I am also attaching the snap shots for dehumidifier geometry and governing equations I have programmed. Could you or @Torsten please help me to trouble shoot the error in m file? Thanks a lot for your help.

% Author: Jubair Shamim, Oak Ridge National Laboratory 
% Date created: Mar 24, 2025
% Date of last modification: TBA

%Model Assumptions
%----------------------------------------------------
% 1.    Heat capacities of SS frame and mesh are neglected 
% 2.    Area of air flow channel is assumed to be same width and length of
%       desiccant bed 


clear all
close all
clc

% Discretization parameters
L = 0.2; % length of air channel in meter (same as desiccant bed)
step=0.002;
N=L/step; % no. of discretization in X direction
X = linspace(0,L,N); % discretise the domain 
n = numel(X);

% Output locations for temperature and humidity 
N_inlet= 2; 
N_middle= N/2; 
N_outlet=N-1; 

% time interval
ti=0;
tend=600; % end time 
delt=1;
Nt=tend/delt;
tspan = linspace (ti,tend,Nt); % given time span
%-------------------------------------------------------------
% 
% CFL= 0.197*(delt/step); % water flow velocity 0.197 m/s
% display (CFL)


% Define initial conditions
Q_init= 1E-10; % initial water in the desiccant bed in kg; 
W_init=1E-10; %initial humidity in the channel in kg/kg_DA; assume 
Tair_init=300.0; %initial air temp in K
Tdes_init=300.0; %initial desiccant temp in K


% Initial conditions input in ODE
Tair0        =  Tair_init *ones(n,1);        
Tdes0        =  Tdes_init  *ones(n,1);     
W0           =  W_init*ones(n,1);     
Q0           =  Q_init*ones(n,1);   

%DC0      = 0     *ones(n,1);     % Conversion rate
%-------------------------------------------------------------

% Initialize output arrays
Q_total = zeros(1, n); % to store total adsorbed water in the bed
W_dehumidification_capacity = zeros(1, n); % to store dehumidification capacity



% setting up the solver
y0 = [Tair0; Tdes0; W0; Q0];
%  y0 = [Tgas0; Q0];
[T,Y] = ode45 (@(t,y) solution(t,y,X,n),tspan,y0);
Y=Y.';
%-------------------------------------------------------------

% saving results
Tair = Y(1  :  n,:);
Tdes = Y(n+1:2*n,:);
W = Y(2*n+1:3*n,:);
Q  = Y(3*n+1:4*n,:);
%DC  = Y(4*n+1:5*n,:);
%-------------------------------------------------------------

%writematrix
writematrix(W, 'W.csv');
writematrix(Tdes, 'Tdes.csv');
writematrix(Tair, 'Tair.csv');
writematrix(Q, 'Q.csv');
%writematrix(DC, 'DC.csv');

%Experiment data
% Table = readtable('Adsorption.xlsx');
% x1 = Table.Time; %: get the excel column, Header1 (header name)
% y1 = Table.Tw_Exp; 
% y2 = Table.T_gas_Exp; 
% y3 = Table.T_bed_Exp; 

% Table = readtable('For_validation_Kimura.xlsx');
% x1 = Table.time; %: get the excel column, Header1 (header name)
% y1 = Table.Two_exp; 
% y2 = Table.Bed_in_Exp; 
% y3 = Table.Bed_Middle_Exp; 
% y4 =Table.Bed_Out_exp; 
% 

f = figure();
f.WindowState = 'maximized';

% Desiccant temperature at inlet, middle and outlet

subplot(2, 2, 1);
plot (T, Tair(N_inlet,1:Nt), 'g', T, Tair (N_middle,1:Nt), 'red',   T, Tair (N_outlet,1:Nt),'b');
xlabel('time (s)')
ylabel('air temperature (K)')   
legend('near inlet', 'at middle', 'near outlet');

subplot(2, 2, 2);
plot (T, Tdes(N_inlet,1:Nt), 'g', T, Tdes (N_middle,1:Nt), 'red',   T, Tdes (N_outlet,1:Nt),'b');
xlabel('time (s)')
ylabel('desiccant temperature (K)')   
legend('near inlet', 'at middle', 'near outlet');

subplot(2, 2, 3);
plot (T, W (N_inlet,1:Nt), 'g', T, W (N_middle,1:Nt), 'red',   T, W (N_outlet,1:Nt),'b');
xlabel('time (s)')
ylabel('Humidity Distribution (Kg/kg)')   
legend('near inlet', 'at middle', 'near outlet')


subplot(2, 2, 4);
plot(T, Q (N/2, 1:Nt));

title('Loading');
xlabel('Time');
ylabel('Loading at middle [kg/kg]');
%----------------------------------------------------------------

function DyDt = solution(~,y,X,n)

% Define Constants 
% MS Gel used as desiccant 

% Inlet parameters 
%---------------------------------------------------
W_in=0.015; % inlet humidity in kg/kg_dry air 
Tair_in = 300E0; % inlet air temperature in K 
U_in=0.9; % inlet air velocity in m/sec.
Tamb= 300; % ambient temperature 


% Geometric properties 
%------------------------------------------------
L_air=0.2; % length of MFBDD air channel, m 

% % Areas for shinge sheet 
% A_des= 1.00E-5; % Cross-section area of single MFBDD desiccant sheet, m2
% A_air= 4.17E-5; % Cross-section area of single MFBDD air channel, m2 (for bed thickness 1 mm)
% A_sa=2E-3; % Heat Transf. surface area between bed and air for 1 sheets in m^2

% Areas for 5 sheets and 6 air channel 
A_des= 5.00E-5; %(Lily)
A_air= 2.5E-4;  %Estimated
A_sa=1E-2; %Estimated
%A_sa=3.4E-2; %(Lily)
MCp= 7; % Mass*Cp of 5 steel frames in J/K (Lily)


% Heat Transfer Properties for air
%--------------------------------------------------
h_sa=7.5E0; %convective HTC between bed and air in W/m^2 K (hsu)
%h_amb_Area=1.0; % HTC*Area between bed and ambient in W/K (Hsu)- 1 sheet 
h_amb_Area=0.5; % HTC*Area between bed and ambient in W/K (Hsu)- 5 sheet 


%Air properties 
%----------------------------------------------------
R=8.3145; %universal gas constant in J/mol. K
M_w= 1.8015E-2; % Mole Mass of water in kg/mol 
rho_air = 1.225E0; %Density of air in kg/m^3
Cp_air =1.027E3; % speific heat of air in J/kg.K


% Desiccant properties 
%--------------------------------------------
poro= 0.3; % bed porosity 
rho_des= 1.1E3; % density of desiccant in kg/m3
Cp_des= 9.21E2; % specific heat of desiccant, J/kg.K

%m_des= 1.52E-3; % mass of desiccant (MS Gel) in 1 sheets in kg 
%m_des= 7.6E-3; % mass of desiccant (MS Gel) in 5 sheets in kg 


% Adsorption properties of MS Gel 
%---------------------------------------------
H_ads=2.26E6; % enthalpy of adsorption in J/kg 
K_LDF  = 1.5E-3; % % Adsorption rate constant in s-1 (Lily)

% Constants for langmuir-sips isotherm 
a= 0.39; 
b1=0.02; 
b2=2E-11;
q=6;


% Heat Transfer Properties for water (TBD)
%--------------------------------------------------


%End of constants------------------------------------------------------------


% Define matrix

%Change in variables with time
DTairDt = zeros(n,1); 
DTdesDt = zeros(n,1); 
DWDt = zeros(n,1); 
DQDt = zeros(n,1); 


% from saving results
Tair = y(1:n);            
Tdes = y(n+1:2*n);        
W = y(2*n+1:3*n);      
Q  = y(3*n+1:4*n);  

%     next variable  = y(4*n+1:5*n);      % Conversion rate


% Constants for Governing equations 

% Mass balance for humidity 
%----------------------------------
C11= (1-poro)*A_des*rho_des;
C12=(A_air+(poro*A_des))*rho_air; 
C1=C11/C12; % for input in eqn 

%Energy Eqn for air 
%------------------------------------
C21=(A_sa*h_sa);
C22=(A_air+(poro*A_des))*rho_air*Cp_air*L_air; 
C2= C21/C22; % for input in eqn 

% Energy Eqn for Bed/desiccant (neglect heat capacity of dehumidifier and frame, to be accounted in loss term)
%----------------------------------------

C31= ((1-poro)*rho_des*Cp_des*A_des*L_air)+MCp; %MCp is the heat caapcity of steel frame
%C31= ((1-poro)*rho_des*Cp_des*A_des*L_air)+10; 
C3= (H_ads*rho_des*(1-poro)*A_des*L_air)/C31; % for input in eqn, heat generation term 

C4= (h_sa*A_sa)/C31; % for input in eqn, HT term with air 

C5= h_amb_Area/C31; % for input in eqn, Heat loss term 



% Governing equations 

%  Q  
for i = 1:n


% Calculate average W and Tair for the discretized domain
    avg_W = mean(W(:,i)); % Average of W for the spatial domain at time step i
    avg_Tair = mean(Tair(:,i)); % Average of Tair for the spatial domain at time step i
    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85)); 
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w; 
    phi = (P_vap * 100) / P_sat; % relative humidity
    % Langmuir-sips isotherm equation
    %Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    Q_eq = ((a*b1*phi)/1+(b1*phi)) + ((a*b2*(phi^q))/(1+(b2*(phi^n)))); % Langmuir-sips isotherm equation 


    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));   
    % Calculate total mass of water adsorbed in the bed (integrated over the domain)
    Q_total(i) = sum(Q) / n; 
    % Calculate dehumidification capacity
    W_inlet = W(1, i); % Inlet (first column in W_data)
    W_outlet = W(n, i); % Outlet (last column in W_data)
    W_dehumidification_capacity(i) = mean(W_outlet) - mean(W_inlet); 
 
end


%Humidity ---------------------------------------------------------------------------------

W(1) = W_in;   % BC at x = 0

for i = 2:n    % interior and last grid point
  
    DWDt(i) =  (-C1*K_LDF*(Q_eq - Q(i))) - (U_in*((W(i)-W(i-1))/(X(i)-X(i-1)))) ;     
   

end  


%air ---------------------------------------------------------------------------------

Tair(1) = Tair_in;   % BC at x = 0

for i = 2:n       % interior and last grid point
  
    DTairDt(i) =  (C2*(Tdes(i) - Tair(i))) - (U_in*((Tair(i)-Tair(i-1))/(X(i)-X(i-1)))) ;     
   
end
%-------------------------------------------------------------------------------------


% Bed/desiccant

for i = 1:n
 DTdesDt(i) = (C3*K_LDF*(Q_eq - Q(i)))- (C4*(Tdes(i)-Tair(i)))- (C5*(Tdes(i)-Tamb)); 
end



DyDt = [DTairDt; DTdesDt;DWDt; DQDt]; 

end   %end of function 

[JUBAIR]

@Torsten, First, we need to estimate phi, which is relative humity, for this we should use average of W and Tair at each time step over the whole domain discretized in n divisions.

For adsorbed mass, Q_total should be the integration of Q over the domain discretized in n divisions at each time step.

There is no missing equation in my m file (SSLC_NoWater_B2_LilyEqn_Copy.m) as atatched. Just getting some error as you can see below. If you could help to remove the error and make the code working, I can do the rest.

Thanks for your help.

Code:

clear all
close all
clc

% Discretization parameters
L = 0.2; % length of air channel in meter (same as desiccant bed)
step=0.002;
N=L/step; % no. of discretization in X direction
X = linspace(0,L,N); % discretise the domain 
n = numel(X);

% Output locations for temperature and humidity 
N_inlet= 2; 
N_middle= N/2; 
N_outlet=N-1; 

% time interval
ti=0;
tend=600; % end time 
delt=1;
Nt=tend/delt;
tspan = linspace (ti,tend,Nt); % given time span
%-------------------------------------------------------------
% 
% CFL= 0.197*(delt/step); % water flow velocity 0.197 m/s
% display (CFL)


% Define initial conditions
Q_init= 1E-10; % initial water in the desiccant bed in kg; 
W_init=1E-10; %initial humidity in the channel in kg/kg_DA; assume 
Tair_init=300.0; %initial air temp in K
Tdes_init=300.0; %initial desiccant temp in K


% Initial conditions input in ODE
Tair0        =  Tair_init *ones(n,1);        
Tdes0        =  Tdes_init  *ones(n,1);     
W0           =  W_init*ones(n,1);     
Q0           =  Q_init*ones(n,1);   

%DC0      = 0     *ones(n,1);     % Conversion rate
%-------------------------------------------------------------

% setting up the solver
y0 = [Tair0; Tdes0; W0; Q0];
%  y0 = [Tgas0; Q0];
[T,Y] = ode45 (@(t,y) solution(t,y,X,n),tspan,y0);
Y=Y.';
%-------------------------------------------------------------

% saving results
Tair = Y(1  :  n,:);
Tdes = Y(n+1:2*n,:);
W = Y(2*n+1:3*n,:);
Q  = Y(3*n+1:4*n,:);
%DC  = Y(4*n+1:5*n,:);
%-------------------------------------------------------------

%writematrix
writematrix(W, 'W.csv');
writematrix(Tdes, 'Tdes.csv');
writematrix(Tair, 'Tair.csv');
writematrix(Q, 'Q.csv');
%writematrix(DC, 'DC.csv');

%Experiment data
% Table = readtable('Adsorption.xlsx');
% x1 = Table.Time; %: get the excel column, Header1 (header name)
% y1 = Table.Tw_Exp; 
% y2 = Table.T_gas_Exp; 
% y3 = Table.T_bed_Exp; 

% Table = readtable('For_validation_Kimura.xlsx');
% x1 = Table.time; %: get the excel column, Header1 (header name)
% y1 = Table.Two_exp; 
% y2 = Table.Bed_in_Exp; 
% y3 = Table.Bed_Middle_Exp; 
% y4 =Table.Bed_Out_exp; 
% 

f = figure();
f.WindowState = 'maximized';

% Desiccant temperature at inlet, middle and outlet

subplot(2, 2, 1);
plot (T, Tair(N_inlet,1:Nt), 'g', T, Tair (N_middle,1:Nt), 'red',   T, Tair (N_outlet,1:Nt),'b');
xlabel('time (s)')
ylabel('air temperature (K)')   
legend('near inlet', 'at middle', 'near outlet');

subplot(2, 2, 2);
plot (T, Tdes(N_inlet,1:Nt), 'g', T, Tdes (N_middle,1:Nt), 'red',   T, Tdes (N_outlet,1:Nt),'b');
xlabel('time (s)')
ylabel('desiccant temperature (K)')   
legend('near inlet', 'at middle', 'near outlet');

subplot(2, 2, 3);
plot (T, W (N_inlet,1:Nt), 'g', T, W (N_middle,1:Nt), 'red',   T, W (N_outlet,1:Nt),'b');
xlabel('time (s)')
ylabel('Humidity Distribution (Kg/kg)')   
legend('near inlet', 'at middle', 'near outlet')


subplot(2, 2, 4);
plot(T, Q (N/2, 1:Nt));

title('Loading');
xlabel('Time');
ylabel('Loading at middle [kg/kg]');
%----------------------------------------------------------------

function DyDt = solution(~,y,X,n)

% Define Constants 
% MS Gel used as desiccant 

% Inlet parameters 
%---------------------------------------------------
W_in=0.015; % inlet humidity in kg/kg_dry air 
Tair_in = 300E0; % inlet air temperature in K 
U_in=0.9; % inlet air velocity in m/sec.
Tamb= 300; % ambient temperature 


% Geometric properties 
%------------------------------------------------
L_air=0.2; % length of MFBDD air channel, m 

% % Areas for shinge sheet 
% A_des= 1.00E-5; % Cross-section area of single MFBDD desiccant sheet, m2
% A_air= 4.17E-5; % Cross-section area of single MFBDD air channel, m2 (for bed thickness 1 mm)
% A_sa=2E-3; % Heat Transf. surface area between bed and air for 1 sheets in m^2

% Areas for 5 sheets and 6 air channel 
A_des= 5.00E-5; %(Lily)
A_air= 2.5E-4;  %Estimated
A_sa=1E-2; %Estimated
%A_sa=3.4E-2; %(Lily)
MCp= 7; % Mass*Cp of 5 steel frames in J/K (Lily)


% Heat Transfer Properties for air
%--------------------------------------------------
h_sa=7.5E0; %convective HTC between bed and air in W/m^2 K (hsu)
%h_amb_Area=1.0; % HTC*Area between bed and ambient in W/K (Hsu)- 1 sheet 
h_amb_Area=0.5; % HTC*Area between bed and ambient in W/K (Hsu)- 5 sheet 


%Air properties 
%----------------------------------------------------
R=8.3145; %universal gas constant in J/mol. K
M_w= 1.8015E-2; % Mole Mass of water in kg/mol 
rho_air = 1.225E0; %Density of air in kg/m^3
Cp_air =1.027E3; % speific heat of air in J/kg.K


% Desiccant properties 
%--------------------------------------------
poro= 0.3; % bed porosity 
rho_des= 1.1E3; % density of desiccant in kg/m3
Cp_des= 9.21E2; % specific heat of desiccant, J/kg.K

%m_des= 1.52E-3; % mass of desiccant (MS Gel) in 1 sheets in kg 
%m_des= 7.6E-3; % mass of desiccant (MS Gel) in 5 sheets in kg 


% Adsorption properties of MS Gel 
%---------------------------------------------
H_ads=2.26E6; % enthalpy of adsorption in J/kg 
K_LDF  = 1.5E-3; % % Adsorption rate constant in s-1 (Lily)

% Constants for langmuir-sips isotherm 
a= 0.39; 
b1=0.02; 
b2=2E-11;
q=6;


% Heat Transfer Properties for water (TBD)
%--------------------------------------------------


%End of constants------------------------------------------------------------


% Define matrix

%Change in variables with time
DTairDt = zeros(n,1); 
DTdesDt = zeros(n,1); 
DWDt = zeros(n,1); 
DQDt = zeros(n,1); 

Q_total = zeros(1, n); % to store total adsorbed water in the bed
W_dehumidification_capacity = zeros(1, n); % to store dehumidification capacity

% from saving results
Tair = y(1:n);            
Tdes = y(n+1:2*n);        
W = y(2*n+1:3*n);      
Q  = y(3*n+1:4*n);  

%Q_total = y(4*n+1:5*n);     
%W_dehumidification_capacity =y(5*n+1:6*n);  

% Constants for Governing equations 

% Mass balance for humidity 
%----------------------------------
C11= (1-poro)*A_des*rho_des;
C12=(A_air+(poro*A_des))*rho_air; 
C1=C11/C12; % for input in eqn 

%Energy Eqn for air 
%------------------------------------
C21=(A_sa*h_sa);
C22=(A_air+(poro*A_des))*rho_air*Cp_air*L_air; 
C2= C21/C22; % for input in eqn 

% Energy Eqn for Bed/desiccant (neglect heat capacity of dehumidifier and frame, to be accounted in loss term)
%----------------------------------------

C31= ((1-poro)*rho_des*Cp_des*A_des*L_air)+MCp; %MCp is the heat caapcity of steel frame
%C31= ((1-poro)*rho_des*Cp_des*A_des*L_air)+10; 
C3= (H_ads*rho_des*(1-poro)*A_des*L_air)/C31; % for input in eqn, heat generation term 

C4= (h_sa*A_sa)/C31; % for input in eqn, HT term with air 

C5= h_amb_Area/C31; % for input in eqn, Heat loss term 



% Governing equations 

%  Q  
for i = 1:n


% Calculate average W and Tair for the discretized domain
    avg_W = mean(W(:,i)); % Average of W for the spatial domain at time step i
    avg_Tair = mean(Tair(:,i)); % Average of Tair for the spatial domain at time step i
    
    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85)); 
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w; 
    phi = (P_vap * 100) / P_sat; % relative humidity
    
    % Langmuir-sips isotherm equation
    %Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    Q_eq = ((a*b1*phi)/1+(b1*phi)) + ((a*b2*(phi^q))/(1+(b2*(phi^n)))); % Langmuir-sips isotherm equation 


    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));   
    
    % Calculate total mass of water adsorbed in the bed (integrated over the domain)
    Q_total(i) = sum(Q) / n; 
    
    
   % Calculate dehumidification capacity
    W_inlet = W(1, i); % Inlet (first column in W_data)
    W_outlet = W(n, i); % Outlet (last column in W_data)
    W_dehumidification_capacity(i) = mean(W_outlet) - mean(W_inlet); 
 
end


%Humidity ---------------------------------------------------------------------------------

W(1) = W_in;   % BC at x = 0

for i = 2:n    % interior and last grid point
  
    DWDt(i) =  (-C1*K_LDF*(Q_eq - Q(i))) - (U_in*((W(i)-W(i-1))/(X(i)-X(i-1)))) ;     
   

end  


%air ---------------------------------------------------------------------------------

Tair(1) = Tair_in;   % BC at x = 0

for i = 2:n       % interior and last grid point
  
    DTairDt(i) =  (C2*(Tdes(i) - Tair(i))) - (U_in*((Tair(i)-Tair(i-1))/(X(i)-X(i-1)))) ;     
   
end
%-------------------------------------------------------------------------------------


% Bed/desiccant

for i = 1:n
 DTdesDt(i) = (C3*K_LDF*(Q_eq - Q(i)))- (C4*(Tdes(i)-Tair(i)))- (C5*(Tdes(i)-Tamb)); 
end



DyDt = [DTairDt; DTdesDt;DWDt; DQDt]; 

end   %end of function 

[Torsten]

The equations you included are simple ODEs, but you try to solve PDEs in your code.

[JUBAIR]

@Torsten the code was working fine before I integrated the following part. Just need to know how to integrate this part correctly.

% Initialize output arrays
Q_total = zeros(1, n); % to store total adsorbed water in the bed
W_dehumidification_capacity = zeros(1, n); % to store dehumidification capacity

for i = 1:n
    % Calculate average W and Tair for the discretized domain
    avg_W = mean(W_data(:, i)); % Average of W for the spatial domain at time step i
    avg_Tair = mean(Tair_data(:, i)); % Average of Tair for the spatial domain at time step i

    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85)); 
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w; 
    phi = (P_vap * 100) / P_sat; % relative humidity

    % Langmuir-sips isotherm equation
    Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    
    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));   

    % Calculate total mass of water adsorbed in the bed (integrated over the domain)
    Q_total(i) = sum(Q) / n; 

    % Calculate dehumidification capacity
    W_inlet = W_data(1, i); % Inlet (first column in W_data)
    W_outlet = W_data(n, i); % Outlet (last column in W_data)
    W_dehumidification_capacity(i) = mean(W_outlet) - mean(W_inlet); 
end

[Torsten]

If you compute Q_eq with mean values for W and T_air over the complete reactor, you don't need to compute Q as dependent on the spatial coordinate. You just need to solve one ordinary differential equation dQ/dt = K_LDF*(Q_eq-Q). Thus the loop "for i = 1:n" is superfluous. But I thought that you wanted to consider the spatial dependence of the adsorbed quantity of water as a function of the spatially dependent quantities W and T, didn't you ?

[JUBAIR]

@Torsten: I want to compute Q as dependent on the spatial coordinate. But Q_eq should not depend on spatial coordinate rather it should be computed using mean values of W and T_air. So we need the loop "for i=1:n" for Q(i). Does this answer your question?

[Torsten]

Look at your loop. K_LDF is constant, you want Q_eq the same and independent of i because you want to compute it with average values for T and W. Thus K_LDF*(Q_eq-Q(i)) will be independent of i because you start with a constant value of 0 for Q in all grid points. So given a time T, the Q-profile over the length of the reactor will be a flat horizontal line.

You only get a profile for Q over the length if you prescribe a profile for Q_eq (or a somehow varying value for K_LDF over the length).

Here is the code to test it:

for i = 1:n
% Calculate average W and Tair for the discretized domain
    avg_W = mean(W); % Average of W for the spatial domain at time step i
    avg_Tair = mean(Tair); % Average of Tair for the spatial domain at time step i
%    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85));
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w;
    phi = (P_vap * 100) / P_sat; % relative humidity
%    % Langmuir-sips isotherm equation
    %Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    Q_eq = ((a*b1*phi)/1+(b1*phi)) + ((a*b2*(phi^q))/(1+(b2*(phi^n)))); % Langmuir-sips isotherm equation
%    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));
end

If you now plot the loading as

subplot(2, 2, 4);
plot (T, Q (N_inlet,1:Nt), 'g', T, Q (N_middle,1:Nt), 'red',   T, Q (N_outlet,1:Nt),'b');
title('Loading');
xlabel('Time');
ylabel('Loading at middle [kg/kg]');

you will see that all three curves are equal.

[JUBAIR]

Btw, what should be the expression for W_avg? Should it be something like below?

For i=1:n W_avg =mean (W(i),1); End

This should return the row vector for W_avg at each time step over the n discertized domain, correct?

[JUBAIR]

@Torsten I have updated the code for Q_eq as below.

%  Q  
for i = 1:n

% Calculate average W and Tair for the discretized domain
    avg_W = mean (W(i)); % Average of W for the spatial domain
    avg_Tair = mean(Tair(i)); % Average of Tair for the spatial domain
    
    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85)); 
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w; 
    phi = (P_vap * 100) / P_sat; % relative humidity
    
    % Langmuir-sips isotherm equation
    %Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    Q_eq = ((a*b1*phi)/1+(b1*phi)) + ((a*b2*(phi^q))/(1+(b2*(phi^n)))); % Langmuir-sips isotherm equation 


    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));   
    
    % Calculate total mass of water adsorbed in the bed (integrated over the domain)
    Q_bed(i) = sum(Q(i))/n; 
    
    
   % Calculate dehumidification capacity
    W_inlet = W(2); % Inlet (first column in W_data)
    W_outlet = W(n); % Outlet (last column in W_data)
    DC(i) = W_inlet - W_outlet; 
 
end

Now I get the results like this:

But Now I facing trouble to save the Q_bed(i) and DC(i). Could you please help me to update the code so that these variables can be saved and plotted? attached is the updated m file.

[Torsten]

Replace the line

function DyDt = solution(~,y,X,n)

by

function [DyDt,Q_bed,DC] = solution(~,y,X,n)

and call the function with the solutions Y(:,i) after the simulation has finished:

[T,Y] = ode45 (@(t,y) solution(t,y,X,n),tspan,y0);
for i = 1:numel(T)
    [~,Q_bed,DC] = solution(T(i),Y(:,i),X,n);
    Q_bed_matrix(i,:) = Q_bed;
    DC_matrix(i,:) = DC;
end

[JUBAIR]

@Marie Anna NOVIELLO and @Torsten Thanks a lot to both of you for your great assistance. The code is working fine now. Here is the result using Avg_W and Avg_Tair to estimate phi . This was the last piece of puzzle to solve. With that I will close this case. Have a great day.

[JUBAIR]

@Torsten I have one more question. I want to run my code first 300 sec for adsorption and next 300 sec for desorption. How to set different values for adsorption and desorption for the following variables in my code?

K_LDF, h_sa, W_in; % inlet humidity in kg/kg_dry air Tair_in; % inlet air temperature in K U_in % inlet air velocity in m/sec.

I also need to change the signs of source terms in the governing equations of mass and energy during adsorption and desorption cycles.

I shall highly appreciate it you can help me to write the code for this part. Thanks 🙏

[Torsten]

Write two functions "solution_adsorption.m" and "solution_desorption.m".

Call ode45 for 300 s with function file "solution_adsorption.m".

Save the results.

Use the results of the last output time of the adsorption phase to define the initial values of the desorption phase.

Call ode45 for 300 s with function file "solution_adsorption.m".

Save the results.

[JUBAIR]

@Torsten this is one way to do it using two separate files. But I want to run it for multiple cycles of adsorption and desorption and want to get a curve for Q (i.e., mass adsorbed in the bed) like shown in the picture below. So, doing it in one single file is convenient. I need a program to alter the values of some variables for sorption and desorption and signs of sources terms in the same file for different time intervals. if you could help I would highly appreciate.

[Torsten]

You will always have to stop and restart integration when switching between the two modes of operation. So before you start including multiple if statements in one file that would cover both adsorption and desorption, the best way is to use two files.

[JUBAIR]

@ thanks for this feedback

[JUBAIR]

@Torsten, I am facing another simple problem in this code. Tried for hours today but couldn't solve. I want to save the values of "phi" at each time step in an excel file. It should be easy but I wasn't successful. Could you please provide the code to save phi in excel file? I shall remain grateful to you.

[Torsten]

Save phi in an array Phi(i) and proceed as suggested for Qbed and DC:

[T,Y] = ode45 (@(t,y) solution(t,y,X,n),tspan,y0);
for i = 1:numel(T)
    [~,Q_bed,DC,Phi] = solution(T(i),Y(:,i),X,n);
    Q_bed_matrix(i,:) = Q_bed;
    DC_matrix(i,:) = DC;
    Phi_matrix(i,:) = Phi;
end

function [DyDt,Q_bed,DC,Phi] = solution(~,y,X,n)
...
for i = 1:n
% Calculate average W and Tair for the discretized domain
    avg_W = mean (W(i)); % Average of W for the spatial domain
    avg_Tair = mean(Tair(i)); % Average of Tair for the spatial domain
    
    % Calculate P_vap and P_sat using average values
    P_sat = 610.78 * exp((17.27 * (avg_Tair - 273.15)) / (avg_Tair - 35.85)); 
    P_vap = (avg_W * R * rho_air * avg_Tair) / M_w; 
    phi = (P_vap * 100) / P_sat; % relative humidity

    % Langmuir-sips isotherm equation
    %Q_eq = ((a * b1 * phi) / (1 + (b1 * phi))) + ((a * b2 * (phi^q)) / (1 + (b2 * (phi^n)))));
    Q_eq = ((a*b1*phi)/1+(b1*phi)) + ((a*b2*(phi^q))/(1+(b2*(phi^n)))); % Langmuir-sips isotherm equation 
    % LDF calculation
    DQDt(i) = K_LDF * (Q_eq - Q(i));   
    
    % Calculate total mass of water adsorbed in the bed (integrated over the domain)
    Q_bed(i) = sum(Q(i))/n; 
    
    
   % Calculate dehumidification capacity
    W_inlet = W(2); % Inlet (first column in W_data)
    W_outlet = W(n); % Outlet (last column in W_data)
    DC(i) = W_inlet - W_outlet; 

    Phi(i) = phi;
end
...
end