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 view (last 30 days)
Show older comments
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 Comments
Aditya
on 17 Nov 2023
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.
Answers (0)
See Also
Categories
Find more on Thermal Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!