Partial Differential Equation Solver with Dynamic Boundary conditions

P0 = 101325; % Atmospheric pressure (Pa)
C1 = 10;
C2 = 62.29;
C3 = 62.29;
C4 = 0.557;
C5 = 2391.3;
C6 = 0.609;

omega_i = 0.008; % Inlet humidity ratio (kg/kg dry air)
Tg_i = 303.15; % Inlet air temperature (K)

% Initial conditions
Ts_init = Tg_i; 
omega_s_init = 0.008; %  solid humidity ratio
omega_o_init = omega_i;
Tg_o_init = Tg_i;

% Sim parameters
N = 180; % Angular resolution
theta = linspace(0, 180, N); % Angular positions in degrees
dt = 1 / (15.5 * 60) / N; % Time step for rotation
tspan = [0, dt * N]; % Simulation time

% ODE system  to access constants
function dYdt = desiccant_wheel_ode(t, Y)
    % Unpack variables
    omega_o = Y(1);
    Ts = Y(2);
    Tg_o = Y(3);
    omega_s = Y(4);

    % Saturation pressure and relative humidity
    Ps = exp(5294 / Ts - 6); % Saturation pressure (Pa)
    phi = (omega_s / 0.622) * P0 / ((0.622 + omega_s) * Ps);

    % Empirical functions
    s1 = (0.24 / 1.5) * (phi^0.333) * (0.622 * P0) / ((0.622 + omega_s)^2 * Ps);
    s2 = (0.24 / 1.5) * (phi^0.333) * (5294 * phi) / (Ts^2);

    % Differential equations
    d_omega_dt = C1 * (omega_i - omega_o) + C2 * (omega_s - omega_o);
    d_Ts_dt = C4 * C5 * (omega_o - omega_s) + C6 * (Tg_o - Ts);
    d_Tg_dt = C1 * (Tg_i - Tg_o) + C3 * (Ts - Tg_o);
    d_omega_s_dt = (s2 / s1) * d_Ts_dt + (C4 / s1) * (omega_o - omega_s);

    % Output derivatives
    dYdt = [d_omega_dt; d_Ts_dt; d_Tg_dt; d_omega_s_dt];
end

% ODE ode45
Y0 = [omega_o_init, Ts_init, Tg_o_init, omega_s_init]; % Initial conditions
[t, Y] = ode45(@desiccant_wheel_ode, tspan, Y0);

% Extract results
omega_o = Y(:, 1) * 1000; % Convert to g/kg dry air
Ts = Y(:, 2) - 273.15; % temp to °C

% Plot data
figure;
plot(theta, omega_o, 'k', 'LineWidth', 2);
xlabel('Degree');
ylabel('Humidity [g/kg dry air]');
title('Adsorption-side Outlet Humidity Ratio');
grid on;

[Torsten]

Can you include the mathematical description of your model ? I don't like to extract it from your code and I don't want to join ResearchGate.

[Emmanuel]

Thanks for responding. I have attached a pdf version of the paper. check at the top

[Torsten]

And could you list the equation numbers in the article that you try to solve ?

[Emmanuel]

I have updated it code but somehow my explanation has disappeared. each time i try to fix something from the code, another simply strings-up. so i had to reduce it to just a single speed 15.5 Rph but it looks like my differential equation setup is not done properly or i am using the wrong solver

[Emmanuel]

the entire setup is based on equation 29,36,32,34 which is was i am trying to solve

[Torsten]

As far as I understand the article, the four equations given are representative for a single cell within the entire tube. Thus you will have to repeat the equations n times (with exchange terms between neighbouring cells) to reproduce the PDE model.

The equations given in the article are already a discretized version of the underlying partial differential equations.

[Emmanuel]

that is exactly what i have been trying for the last 24hrs and sadly no luck

[Torsten]

I think the main problem are the 4*n equations themselves, not their implementation in MATLAB. Or could you derive the system you have to solve from the article ? I was not even able to understand the underlying physical system. What is the spatial variable in the direction of which the system is discretized ?

[Emmanuel]

The length of the rotor L (axial direction)

[Torsten]

And the angular position is computed from the axial position via the rotational speed ? Or how does the "degree" variable come into play in the equations ?

[Emmanuel]

yes

[Torsten]

So 180 degrees does not always correspond to axial position L ? Since if rotational speed is higher/lower, angular position of 180 degrees should be reached before/after length L.

[Emmanuel]

So the angular position is computed from the axial position via the rotational speed. theta = omega*t which can be converted to degree. time = the axial position/ axial air velocity. Atleast to my understanding this is how it should be

[Emmanuel]

correct. 180 doesn't always corresponds to L because the angular position and axial position depends on the rotational speed of the wheel. so at high speed, the angular position 180° is technically at some point where the axial position is < L and vice versa for slower speed.

[Torsten]

But at L, the air leaves the apparatus - thus rotational speeds that give 180° after L is reached wouldn't make sense. Further, I don't think that graphs for different contact times can be compared.

I'm in doubt that this is really the correct interpretation of the results.

[Emmanuel]

This is correct cos after L is reached, any further rotational speed will just be meaningless to the system. So take the case of Fig5B, the speed is optimized so when the angle is at 180, the axial position is exactly = L.

[Torsten]

But then, why simulations are necessary ? If L is given, rotational speed can immediately be adjusted such that at 180°, the outlet L is reached.

[Emmanuel]

I don't understand u fully. Whilst it seems straightforward both the humidity and temp profiles are from dynamic mass and heat transfer processes which depends on boundary/inital conditions and other parameters. They are also time dependent

[Torsten]

If you neglect the diffusive/conductive terms in Table 1 listed in this document

https://www.irbnet.de/daten/iconda/CIB8396.pdf

you get the relevant equations before discretization in axial direction. These are the basis for a solution outside SIMULINK.

[Emmanuel]

even more confusing to me. Can u pls provide an example code to this. It really will be so help

[Torsten]

Look up "method-of-lines" as a method to solve partial differential equations.

The method aims at replacing the spatial derivative (d/dx) in the equations by finite difference quotients and leaving the temporal derivatives in their continuous form. The resulting system of ordinary differential equations in the grid points can then be solved by ode45, e.g..

In principle, in the original article, this discretization is done. The equations written therein are the equations in a single grid point in x-direction.

I found the new article could clarify many open questions (at least for me) - especially the question of how the apparatus works.

[Torsten]

Since the results from the first article are plotted over the angle, I suspect that the spatial coordinate according to which the discretization is done there is the angle theta, not the axial length x. But I'm not going to dive into it more deeply.

[Emmanuel]

Well, thanks for your import. Much appreciated though i am still trying to figure work. I think the problem is with the discreization cos i found similar paper which models just a normal Desiccant unit with as different speeds(Table 2). Still, everything get complicated pretty fast.

[Torsten]

Here is a publication where both spatial directions (length and angles) are taken into account.

[Emmanuel]

I tried going back to basics by using the equations on the paper you provided but to validate experimental results 2 on this paper which is attached here. Cos if i can get this right then i can update my previous code successfully. still troubles here again. Results are not looking good

L = 0.2;          % Wheel length (m)
phi = 0.0349;     % Rotational speed (rad/s) for 20 RPH
ua = 1.5;         % Process air velocity (m/s)
ur = 2.5;         % Regeneration air velocity (m/s)

% Solved/Given DW Properties
rho_g = 1.13;     % Air density (kg/m³)
cp_g = 1007;      % Air specific heat (J/kg·K)
ht = 40.8;        % Heat transfer coefficient (W/m²·K)
hm = 0.037;       % Mass transfer coefficient (m/s)
fv_fs = 2.33e-3;  % Volume to surface ratio (m)
rho_f = 1129;     % Desiccant density (kg/m³)
cp_f = 921;       % Desiccant specific heat (J/kg·K)
kf = 0.078;       % Desiccant thermal conductivity (W/m·K)
Ds = 1.2e-6;      % Surface diffusion coefficient (m²/s)
q = 2897e3;       % Heat of adsorption (J/kg)
delta = 2e-4;     % Thickness of desiccant (m)

% Input conditions (Case 16)
Tap_in = 35.9;    % Process air inlet temp (°C)
Tar_in = 120.4;   % Regeneration air inlet temp (°C)
Yap_in = 0.0067;  % Process air inlet humidity (kg/kg)

% Numerical parameters
Nz = 50;          % Number of axial points
Nt = 1000;        % Number of time steps
dz = L/Nz;
dt = 0.01;        % Time step (s)

% Ini arrays
T_g = zeros(Nz, Nt);
T_f = zeros(Nz, Nt);
Y_g = zeros(Nz, Nt);
Y_f = zeros(Nz, Nt);

% Ini conditions
T_g(:,1) = Tap_in;
T_f(:,1) = Tap_in;
Y_g(:,1) = Yap_in;
Y_f(:,1) = 0.045;

% Time stepping
for t = 1:Nt-1
    for z = 2:Nz-1
        % terms for air
        PHI_T = ht*(fv_fs)*(T_f(z,t) - T_g(z,t));
        PHI_TM = rho_f*hm*(fv_fs)*cp_g*T_g(z,t)*(Y_f(z,t) - Y_g(z,t));
        PHI_M = rho_f*hm*(fv_fs)*(Y_f(z,t) - Y_g(z,t));
        
        % terms for desiccant
        PSI_T = ht*(delta/2)*(T_g(z,t) - T_f(z,t));
        PSI_TM = rho_f*q*hm*(delta/2)*(Y_g(z,t) - Y_f(z,t));
        PSI_M = rho_f*hm*(delta/2)*(Y_g(z,t) - Y_f(z,t));
        
        % Air equations
        dTg_dt = -ua/(rho_g*cp_g) * (T_g(z+1,t) - T_g(z-1,t))/(2*dz) + ...
                 (PHI_T + PHI_TM)/(rho_g*cp_g);
        
        dYg_dt = -ua/rho_g * (Y_g(z+1,t) - Y_g(z-1,t))/(2*dz) + ...
                 PHI_M/rho_g;
        
        % Desiccant equations
        dTf_dt = kf/(rho_f*cp_f) * (T_f(z+1,t) - 2*T_f(z,t) + T_f(z-1,t))/dz^2 + ...
                 (PSI_T + PSI_TM)/(rho_f*cp_f);
        
        dYf_dt = Ds * (Y_f(z+1,t) - 2*Y_f(z,t) + Y_f(z-1,t))/dz^2 + ...
                 PSI_M/rho_f;
        
        % Update values
        T_g(z,t+1) = T_g(z,t) + dt*dTg_dt;
        Y_g(z,t+1) = Y_g(z,t) + dt*dYg_dt;
        T_f(z,t+1) = T_f(z,t) + dt*dTf_dt;
        Y_f(z,t+1) = Y_f(z,t) + dt*dYf_dt;
    end
    
    % Boundary conditions
    T_g(1,t+1) = Tap_in;
    Y_g(1,t+1) = Yap_in;
    T_g(Nz,t+1) = T_g(Nz-1,t+1);
    Y_g(Nz,t+1) = Y_g(Nz-1,t+1);
    
    % wheel rotation
    if mod(t*dt, 2*pi/phi) < dt
        T_g(:,t+1) = flipud(T_g(:,t+1));
        Y_g(:,t+1) = flipud(Y_g(:,t+1));
    end
end

% output conditions
rev_steps = min(round(2*pi/(phi*dt)), Nt);
Tap_out = mean(T_g(Nz,max(1,end-rev_steps+1):end));
Yap_out = mean(Y_g(Nz,max(1,end-rev_steps+1):end));

fprintf('Tap_out = %.2f °C\n', Tap_out);
fprintf('Yap_out = %.4f kg/kg\n', Yap_out);

[Torsten]

Use

        % Air equations
        dTg_dt = -ua/(rho_g*cp_g) * (T_g(z,t) - T_g(z-1,t))/(dz) + ...
                 (PHI_T + PHI_TM)/(rho_g*cp_g);
        
        dYg_dt = -ua/rho_g * (Y_g(z,t) - Y_g(z-1,t))/(dz) + ...
                 PHI_M/rho_g;
        

instead of

                % Air equations
        dTg_dt = -ua/(rho_g*cp_g) * (T_g(z+1,t) - T_g(z-1,t))/(2*dz) + ...
                 (PHI_T + PHI_TM)/(rho_g*cp_g);
        
        dYg_dt = -ua/rho_g * (Y_g(z+1,t) - Y_g(z-1,t))/(2*dz) + ...
                 PHI_M/rho_g;
        

and

        dt = 0.001;        % Time step (s)
        

instead of

        dt = 0.01;        % Time step (s)

[Emmanuel]

Thanks for the corrections

[Emmanuel]

The dynamic process is still not captured. Normally, Yap_out is suppose to decrease, but it is increasing quiet significant. Same with Tap_out which is no were near the expected outcome.

Tap_out = 35.96 °C

Yap_out = 0.0071 kg/kg