% This program calculates distillation lines for specific composition points for ternary mixtures.
% Isobaric bubble point calculations  
% VAPOR-LIQUID EQUILIBRIUM for nonideal multicomponent mixtures 
%                                                      
%      y_i P = x_i Psat_i(T) gamma_i(T,x)              
%                                                      
% 1. Specify P and x              
% 2. Guess T                                                   
% 3. Calculate Psat_i, i=1,...,nc from Antoine equation      
% 4. Calculate gamma_i, i=1,...,nc from NRTL equation 
%    (or another liquid activity coefficient model) 
% 5. Calculate y_i from VLE equation above  
% 6. If sum(y_i) is sufficiently close to unity, then stop.      
%    Otherwise, adjust T and go to step 3.

% Ternary mixture of comp1, comp2 and comp3
% Pure components boiling points, Tbp=[80.1 100 101.3]=[353.25  373.15  374.45]

global X P NC R W A alpha 

R=1.98721;   % universal gas constant [cal/molK]

NC=3;   % number of components
comp1='Benzene';
comp2='Water';
comp3='1,4-Dioxane';

P=1;    % pressure in atm

% Data from the DEECHEMA data series (Gmehling and Onken, 1977)
% Antoine coefficients where Psat_i [mmHg] and T [K]
A=[6.87987  1196.760  219.161     
   8.07131  1730.630  233.426  
   7.43155  1554.679  240.337];

% NRTL parameters (small parameters, Wij) 
W=[   0        1522.2402   -293.8487
   1821.8322     0         1551.5163
    434.1172   1097.8251      0     ];

% Alpha parameter
alpha=[0        0.3547      0.3022
       0.3547      0        0.5457      
       0.3022   0.5457   0     ];


% Creates the ternary diagram
X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 
hold on
axis off
axis([0 2*sqrt(3)/3 0 1])
axis equal
X = [0 sqrt(3)/3];
Y = [0 1];
plot(X,Y,'k') 
X = [sqrt(3)/3 2*sqrt(3)/3]; 
Y = [1 0];  
h=gca;
text(sqrt(3)/3-0.095,1.07,'Benzene [s]'); h2=get(h,'Children');set(h2(1),'FontWeight','bold'); 
text(sqrt(3)/3-0.055,1.035,'80.1 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(2*sqrt(3)/3-0.08,-0.04,'Water [sn]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.06,-0.08,'100.0 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(-0.14,-0.04,'1,4-Dioxane [sn]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(-0.060,-0.08,'101.3 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

text(2*sqrt(3)/3-0.37,0.71,'het.az [un]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.33,0.66,'69.0 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(2*sqrt(3)/3-0.64,-0.04,'hom.az [s]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.61,-0.08,'86.6 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

text(2*sqrt(3)/3-1.30,0.85,'Serafimov^,s class 2.0-2b');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-1.40,0.70,'____ distillation boundary ');h2=get(h,'Children');set(h2(1),'FontWeight','bold', 'color','b');
text(2*sqrt(3)/3-1.40,0.60,'____ batch distillation boundary');h2=get(h,'Children');set(h2(1),'FontWeight','bold','color','r');
text(2*sqrt(3)/3-1.40,0.50,'- - - - distillation lines');h2=get(h,'Children');set(h2(1),'FontWeight','bold', 'color','g');

plot(X,Y,'k') 

% --------------------------------------------
% Sketch the batch distillation boundary (running from the heteroazeotrope to the dioxane vertex)
% --------------------------------------------
XL1=[0.6642;0];
XL2=[0.3358;0];
tf1=XL1;
tf3=1-XL1-XL2;
tf2=2*sqrt(3)/3 - sqrt(3)/3*tf1 - 2*sqrt(3)/3*tf3;
plot(tf2,tf1,'r-')


% Initial temperatures to start the bubble point calculations (boiling point of water)
T_K=373.15;

% Prepare xy equilibrium vector 
txy1=ones(2,101)';
txy2=ones(2,101)';

Y=[0.40 0.002 0.597]';  % initial value  
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'g-.')


Y=[0.015 0.001 0.984]';  
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'g-.')


Y=[0.0000001 0.0000001 0.9999998]'; 
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'b-')   % Distillation boundary running from the heteroazeotrope to the homoazeotrope 


%------------------------------------------------------------------
Y=[0.005 0.97 0.025]';  
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'g-.')