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
4 次查看(过去 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 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!