I am trying to write the MATLAB code of the question ofresearch paper Drive link of paper-https://drive.google.com/file/d/12N5ugCymHKBubfqOlM2LU3DVBxdrCoUR/view?usp=drive_link
1 次查看(过去 30 天)
显示 更早的评论
function batch_distillation()
% Binary Data
ki = [1.51; 1.13; 0.95; 0.91];
alphaij = [5.91, 2.51, 1.29, 1.00;
1/5.91, 1, 1/2.51, 1/1.29;
1/2.51, 1/2.51, 1, 1/1.29;
1/1.29, 1/1.29, 1/1.29, 1];
% Given data
F_feed = 400; % kmol/day
F_i = [100; 100; 100; 100]; % kmol of each component in the feed
x_i = [0.97; 0.97; 0.97; 0.97]; % molar fractions in the feed
K = ki .* x_i; % Equilibrium constants
alpha = alphaij .* (x_i ./ x_i'); % Relative volatilities
% Time parameters
t_span = [0 10]; % Set the time span for batch distillation (adjust as needed)
t_steps = 100; % Number of time steps
% Initial conditions
L0 = 0; % Initial liquid holdup
V0 = 0; % Initial vapor holdup
z_L0 = x_i; % Initial liquid composition (same as feed composition)
z_V0 = x_i; % Initial vapor composition (same as feed composition)
% Solve the batch distillation using ODE solver
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[t, y] = ode45(@(t,y) batch_distillation_ode(t, y, F_feed, ki, alphaij, x_i, t_span), t_span, [L0; V0; z_L0; z_V0], options);
% Extract results
L = y(:, 1);
V = y(:, 2);
z_L = y(:, 3:6);
z_V = y(:, 7:10);
% Plot results
figure;
plot(t, L, 'b', 'LineWidth', 2);
hold on;
plot(t, V, 'r', 'LineWidth', 2);
xlabel('Time (hours)');
ylabel('Moles (kmol)');
legend('Liquid holdup (L)', 'Vapor holdup (V)');
title('Moles vs Time');
figure;
plot(t, z_L, 'LineWidth', 2);
hold on;
plot(t, z_V, 'LineWidth', 2);
xlabel('Time (hours)');
ylabel('Molar Fraction');
legend('Benzene (L)', 'Toluene (L)', 'Ethylbenzene (L)', 'o-xylene (L)',...
'Benzene (V)', 'Toluene (V)', 'Ethylbenzene (V)', 'o-xylene (V)');
title('Molar Fractions vs Time');
end
function dydt = batch_distillation_ode(t, y, F_feed, ki, alphaij, x_i, t_span)
% Unpack variables
L = y(1);
V = y(2);
z_L = y(3:6);
z_V = y(7:10);
% Calculate liquid and vapor flow rates
F_L = L / t_span(end);
F_V = V / t_span(end);
% Calculate liquid and vapor compositions
x_L = z_L ./ sum(z_L);
x_V = z_V ./ sum(z_V);
% Calculate liquid and vapor phase rates
r_L = F_L .* x_L;
r_V = F_V .* x_V;
% Calculate material balances
dLdt = F_feed - F_L;
dVdt = F_L - F_V;
dz_Ldt = (r_L - z_L) ./ L;
dz_Vdt = (r_V - z_V) ./ V;
% Pack the derivatives
dydt = [dLdt; dVdt; dz_Ldt; dz_Vdt];
end
3 个评论
Aditya
2023-11-17
It seems that due to your initialization all the values are coming to be INF for z_L and z_V, please do check your initialization. if the issue still persists please do share the paper that you are refering to.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Thermal Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!