function Y=NRTL(xA,xB,T,NT)
%*******************************************************************************
% Program written by Stathis Skouras, September 2018
%*******************************************************************************
% Program for ternary VLE by using NRTL equation
%*******************************************************************************
% Ternary system
% (1): Trioxane (light)
% (2): Water (intermediate)
% (3): Mmim Dmp (heavy)
%*******************************************************************************
% Pressure is 1 atm = 1.013 bar = 760 mmHg
% Composition vectors xA, xB, are line vectors of NT columns
%*******************************************************************************
% Gas constant
R=1.98721; % cal/(mol*K)
%*******************************************************************************
% Data that have to be changed for different mixtures
% Pure component constants for Antoine equation (Gmehling et al, 1977-1990)
A1=8.08097; B1=1582.271; C1=239.726; % (T-range: -15C - 84C, Component 1)
A2=8.11220; B2=1592.864; C2=226.184; % (T-range: -20C - 93C, Component 2)
A3=8.37895; B3=1788.020; C3=227.438; % (T-range: -15C - 98C, Component 3)
% Interaction parameters Aij for NRTL equation (Gmehling et al. 1977-1990)
A11=0; % Pure component 1 (cal/mol)
A22=0; % Pure component 2 (cal/mol)
A33=0; % Pure component 3 (cal/mol)
A12=67.2902; % Binary mixture 1-2 (cal/mol)
A21=-70.5092; % Binary mixture 1-2 (cal/mol)
A13=144.4797; % Binary mixture 1-3 (cal/mol)
A31=-12.7427; % Binary mixture 1-3 (cal/mol)
A23=-2.5594; % Binary mixture 2-3 (cal/mol)
A32=56.2391; % Binary mixture 2-3 (cal/mol)
% Alpha interaction parameter
alpha11=0;
alpha22=0;
alpha33=0;
alpha12=0.3009; % Binary mixture 1-2 (cal/mol)
alpha21=alpha12; % Binary mixture 2-1 (cal/mol)
alpha13=0.3067; % Binary mixture 1-3 (cal/mol)
alpha31=alpha13; % Binary mixture 3-1 (cal/mol)
alpha23=0.3007 ; % Binary mixture 2-3 (cal/mol)
alpha32=alpha23; % Binary mixture 3-2 (cal/mol)
% End of data that have to be changed for different mixtures
%*******************************************************************************
%Vapor pressures by Antoine equation
P1s=10.^(A1-B1./(T'-273.15+C1)); % Component 1 (mmHg)
P2s=10.^(A2-B2./(T'-273.15+C2)); % Component 2 (mmHg)
P3s=10.^(A3-B3./(T'-273.15+C3)); % Component 3 (mmHg)
% Binary interaction parameters 'tafji' for NRTL equation
taf11=A11./(R*T'); % Pure component 1
taf22=A22./(R*T'); % Pure component 2
taf33=A33./(R*T'); % Pure component 3
taf12=A12./(R*T'); % Binary mixture 1-2
taf21=A21./(R*T'); % Binary mixture 1-2
taf13=A13./(R*T'); % Binary mixture 1-3
taf31=A31./(R*T'); % Binary mixture 1-3
taf23=A23./(R*T'); % Binary mixture 2-3
taf32=A32./(R*T'); % Binary mixture 2-3
% Binary interaction parameters 'Gji' for NRTL equation
G11=exp(-alpha11*taf11); % Pure component 1
G22=exp(-alpha22*taf22); % Pure component 2
G33=exp(-alpha33*taf33); % Pure component 3
G12=exp(-alpha12*taf12); % Binary mixture 1-2
G21=exp(-alpha21*taf21); % Binary mixture 1-2
G13=exp(-alpha13*taf13); % Binary mixture 1-3
G31=exp(-alpha31*taf31); % Binary mixture 1-3
G23=exp(-alpha23*taf23); % Binary mixture 2-3
G32=exp(-alpha32*taf32); % Binary mixture 2-3
i=1:NT;
% Logarithm of activity coefficients 'LNGAMMA'
LNGAMMA1(i)=(taf11.*G11.*xA(i)+taf21.*G21.*xB(i)+taf31.*G31.*(1-xA(i)-xB(i)))./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i)))...
+ xA(i).*G11./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))).*(taf11 - (xA(i).*taf11.*G11+xB(i).*taf21.*G21+(1-xA(i)-xB(i)).*taf31.*G31)./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))))...
+ xB(i).*G12./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))).*(taf12 - (xA(i).*taf12.*G12+xB(i).*taf22.*G22+(1-xA(i)-xB(i)).*taf32.*G32)./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))))...
+ (1-xA(i)-xB(i)).*G13./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))).*(taf13 - (xA(i).*taf13.*G13+xB(i).*taf23.*G23+(1-xA(i)-xB(i)).*taf33.*G33)./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))));
LNGAMMA2(i)=(taf12.*G12.*xA(i)+taf22.*G22.*xB(i)+taf32.*G32.*(1-xA(i)-xB(i)))./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i)))...
+ xA(i).*G21./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))).*(taf21 - (xA(i).*taf11.*G11+xB(i).*taf21.*G21+(1-xA(i)-xB(i)).*taf31.*G31)./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))))...
+ xB(i).*G22./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))).*(taf22 - (xA(i).*taf12.*G12+xB(i).*taf22.*G22+(1-xA(i)-xB(i)).*taf32.*G32)./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))))...
+ (1-xA(i)-xB(i)).*G23./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))).*(taf23 - (xA(i).*taf13.*G13+xB(i).*taf23.*G23+(1-xA(i)-xB(i)).*taf33.*G33)./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))));
LNGAMMA3(i)=(taf13.*G13.*xA(i)+taf23.*G23.*xB(i)+taf33.*G33.*(1-xA(i)-xB(i)))./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i)))...
+ xA(i).*G31./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))).*(taf31 - (xA(i).*taf11.*G11+xB(i).*taf21.*G21+(1-xA(i)-xB(i)).*taf31.*G31)./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))))...
+ xB(i).*G32./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))).*(taf32 - (xA(i).*taf12.*G12+xB(i).*taf22.*G22+(1-xA(i)-xB(i)).*taf32.*G32)./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))))...
+ (1-xA(i)-xB(i)).*G33./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))).*(taf33 - (xA(i).*taf13.*G13+xB(i).*taf23.*G23+(1-xA(i)-xB(i)).*taf33.*G33)./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))));
% Activity coefficients 'GAMMA'
GAMMA1(i)=exp(LNGAMMA1(i)); % Component 1
GAMMA2(i)=exp(LNGAMMA2(i)); % Component 2
GAMMA3(i)=exp(LNGAMMA3(i)); % Component 3
% Total pressure in each stage (bar)
DP=1-(xA(i).*GAMMA1(i).*P1s(i) + xB(i).*GAMMA2(i).*P2s(i)+(1-xA(i)-xB(i)).*GAMMA3(i).*P3s(i))/760;
% Vapor Liquid Equilibria (VLE)
yA(i)=xA(i).*GAMMA1(i).*P1s(i)./(xA(i).*GAMMA1(i).*P1s(i) + xB(i).*GAMMA2(i).*P2s(i) + (1-xA(i)-xB(i)).*GAMMA3(i).*P3s(i)); % Component 1
yB(i)=xB(i).*GAMMA2(i).*P2s(i)./(xA(i).*GAMMA1(i).*P1s(i) + xB(i).*GAMMA2(i).*P2s(i) + (1-xA(i)-xB(i)).*GAMMA3(i).*P3s(i)); % Component 2
yC(i)=(1-xA(i)-xB(i)).*GAMMA3(i).*P3s(i)./(xA(i).*GAMMA1(i).*P1s(i) + xB(i).*GAMMA2(i).*P2s(i) + (1-xA(i)-xB(i)).*GAMMA3(i).*P3s(i)); % Component 3
Y=[yA;yB;T';DP];
评论9