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;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; 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; desiccant wheel MATLAB Answers — New Questions
