Antoine Class for Vapor-Liquid Equilibrium CalculationsThe Antoine class provides a computational framework to help doing simple vapor-liquid equilibrium calculations in Matlab. These notes introduce Antoine's equation, then shows how the Antoine class can be used to do typical vapor-liquid equilibrium calculations.
These are a few functions and variables useful in formatting the output of subsequenct calculations. clear all;clc;% Utility functions for formatting computational output.sep = '------------------------------------------------------------------';problem = @(s) disp(sprintf('\n%s\nProblem: %s\n%s',sep,s,sep));answer = @(str,val,units) disp(sprintf('%s = %6.3g [%s]',str,val,units));Review: Saturation Pressure of a Pure ComponentAntoine's equation is a representation of the vapor-liquid equilibrium for a pure component. A common form of the equation is where coefficients A,B, and C are determined from experimental data. Values of these coefficients for a large number of chemical species can be found from many sources, including http://webbook.nist.gov/chemistry. Commonly used units are [mmHg] for pressure and [deg C] for temperature. Here we demonstrate use of Antoine's equation to plot the saturation pressure of pure water as a function of temperature. % For water between 60 and 150 degrees CA = 7.96681;B = 1668.21;C = 228;% Temperature rangeT = 60:150;% Antoine's equationlog10P = A - B./(T+C);% Computing PPsat = 10.^log10P;% Construct the desired plotplot(T,Psat);xlabel('Temperature [C]');ylabel('Pressure [mmHg]');title('Saturation Pressure from Antoine Equation');Antoine Class: Object Oriented Vapor/Liquid CalculationsLet's now introduce the Antoine class object. Recent versions of Matlab include methods for objected oriented programming. These methods encapsulate data and functions for in convenient tools for engineering calculations. Objects provide an elegant way to implement engineering calculations once you get used to the concepts and syntax, The Antoine Class DatasetThe Antoine class includes parameters for a small number of compounds. These can be accessed as indicated below. % Load available data into a structure pp = Antoine.data;% Show list of available compoundsproblem('Antoine Class Dataset');disp(p);% Display Antoine data for one member of the datasetproblem('Antoine Data for Ammonia');disp(p.nh3)disp(sep);------------------------------------------------------------------Problem: Antoine Class Dataset------------------------------------------------------------------ h2o: [1x1 Antoine] n2: [1x1 Antoine] o2: [1x1 Antoine] nh3: [1x1 Antoine] meoh: [1x1 Antoine] etoh: [1x1 Antoine] bnz: [1x1 Antoine] tol: [1x1 Antoine] CH4: [1x1 Antoine] C2H6: [1x1 Antoine] C3H8: [1x1 Antoine] nC4: [1x1 Antoine] nC5: [1x1 Antoine] nC6: [1x1 Antoine] nC7: [1x1 Antoine] nC8: [1x1 Antoine]------------------------------------------------------------------Problem: Antoine Data for Ammonia------------------------------------------------------------------ Antoine Properties: Tmin: -83 Tmax: 60 A: 7.3605 B: 926.1320 C: 240.1700 Pmin: 29.3732 Pmax: 1.8843e+04------------------------------------------------------------------Creating Objects for Compounds not in the Antoine Class DatasetBecause only a small number of compounds are included in the dataset included with the Antoine class, users will generally need to create Antoine objects using data gleaned from other sources. The syntax for creating an Antoine object is x = Antoine(Tmin, Tmax, A, B, C)where [Tmin,Tmax] is the temperature range, and A, B, C are the Antoine coefficients. The following codes create two Antoine objects corresponding to benzene and toluene, respectively. Data comes from Table B.4 of Murphy's textbook. bnz = Antoine( 8, 113, 6.90656, 1211.033, 220.79)tol = Antoine( 6, 137, 6.95464, 1344.8, 219.48)bnz = Antoine Properties: Tmin: 8 Tmax: 113 A: 6.9066 B: 1.2110e+03 C: 220.7900 Pmin: 41.0537 Pmax: 1.8986e+03tol = Antoine Properties: Tmin: 6 Tmax: 137 A: 6.9546 B: 1.3448e+03 C: 219.4800 Pmin: 9.7831 Pmax: 1.5212e+03Piecewise specification of Antoine ParametersSometimes the Antoine parameters are given in 'piecewise' fashion over multiple temperature intervals. The Antoine class will handle these cases provided the temperature intervals are (1) contiguous and (2) non-overlapping. An example h2o = Antoine([ 0, 60, 8.10785, 1750.286, 235.0; ... 60, 150, 7.96681, 1668.210, 228.0])h2o.plot;h2o = Antoine Properties: Tmin: [2x1 double] Tmax: [2x1 double] A: [2x1 double] B: [2x1 double] C: [2x1 double] Pmin: [2x1 double] Pmax: [2x1 double]Searching the NIST Chemistry WebbookThe NIST Chemistry Webbook (http://webbook.nist.gov/chemistry/) is an excellent source of thermodynamic data. The Antoine.search method streamlines searching the webbook, and the Antoine.SI creates an Antoine object where parameters are in SI units (K, bar) rather than the more traditional (C, mmHg). % Search for Antoine parameters for propaneAntoine.search('butane');% Create Antoine object for butane given parameters in SI unitsnC4 = Antoine.SI([ 135.42, 212.89, 4.70812, 1200.475, -13.013; ... 212.89, 272.66, 3.85002, 909.650, -36.146; ... 272.66, 425, 4.35576, 1175.581, -2.071]);Psat: Method for Computing Saturation Pressure% The Antoine class includes methods for several standard computations.% This demonstrate calculation of saturation pressure for benzene and% toluene using the |Psat| method.T = 70:110;plot(T, bnz.Psat(T), T, tol.Psat(T));legend('Benzene','Tolune','Location','NW');xlabel('Temperature [deg C]');ylabel('Pressure [mmHg]');title('Saturation Pressure');Tsat: Method for Computing Saturation TemperatureThe Tsat method computes that saturation temperature as a function of pressure. This is a convenient way to compute normal boiling points. % Retrieve Antoine parameters for waterp = Antoine.data;h2o = p.h2o;% Use Tsat method to report boiling point of waterproblem('Normal Boiling Points');answer('Normal Boiling Point of Water',p.h2o.Tsat(760),'deg C');disp(sep);% Use Tsat method to report normal boiling points for all compounds in the% Antoine class data set.disp('Normal Boiling Points of Species in Antoine''s Database');f = fields(p);for i = 1:length(f) disp(sprintf('%6s %7.2f [deg C]',f{i},p.(f{i}).Tsat(760)))enddisp(sep)------------------------------------------------------------------Problem: Normal Boiling Points------------------------------------------------------------------Normal Boiling Point of Water = 100 [deg C]------------------------------------------------------------------Normal Boiling Points of Species in Antoine's Database h2o 100.00 [deg C] n2 -195.42 [deg C] o2 -182.85 [deg C] nh3 -33.43 [deg C] meoh 64.71 [deg C] etoh 78.33 [deg C] bnz 80.03 [deg C] tol 110.63 [deg C] CH4 -161.45 [deg C] C2H6 -88.67 [deg C] C3H8 -42.73 [deg C] nC4 -0.38 [deg C] nC5 35.06 [deg C] nC6 68.73 [deg C] nC7 98.43 [deg C] nC8 125.68 [deg C]------------------------------------------------------------------plot: Plotting Saturation Pressure versus TemperatureThe Antoine class includes a simple plotting function for saturation pressure. p = Antoine.data;clf;hold on;p.n2.plot([],'b');p.o2.plot([],'r');hold off;legend('Nitrogen','Oxygen','Location','NW')gridApplication: Normal Boiling Point of a Pure ComponentThe boiling point of a pure component is the temperature at which the saturation pressure is equal to the ambient pressure. A simple way to find the boiling point is to use Matlab's fzero function to solve Antoine's equation for temperature as a function of pressure. To use fzero, we need a function which evaluates to 0 for the desired temperature and pressure. Here we use an anonymous function for this purpose (if you haven't encountered these before, anonymous functions are useful Matlab technique for representing and manipulating functions that can be represented by a single line of Matlab code). % Antoine's equation in a form f(T,P) = 0f = @(P,T) log10(P) - A + B/(T+C);% The normal boiling point of water is the solution to f(T,760) = 0Tb = fzero(@(T)f(760,T),[60,150]);problem('Normal Boiling Point of Water');answer('Normal boiling point of water (fzero) ',Tb,'deg C');% The |Tsat| method provides a direct method for computing the saturation% temperature as a function of pressure.p = Antoine.data;answer('Normal boiling point of water (Tsat method)',p.h2o.Tsat(760),'deg C');disp(sep);------------------------------------------------------------------Problem: Normal Boiling Point of Water------------------------------------------------------------------Normal boiling point of water (fzero) = 100 [deg C]Normal boiling point of water (Tsat method) = 100 [deg C]------------------------------------------------------------------Application: What is the boiling point of water in Littleton, Colorado?The average barometric pressure in Littleton, Colorado, is 24.63 in of Hg (compared to sea level 29.92 in Hg). What is the boiling point of water at this reduced pressure? % Convert inHg to mmHg (25.4 mm/in)P = 24.63*25.4;% Boiling point of waterproblem('Boiling Point of Water in Littleton, Colorado');answer('Atmospheric Pressure',P,'mm Hg');answer(' Boiling Point', p.h2o.Tsat(P),'deg C');disp(sep);------------------------------------------------------------------Problem: Boiling Point of Water in Littleton, Colorado------------------------------------------------------------------Atmospheric Pressure = 626 [mm Hg] Boiling Point = 94.6 [deg C]------------------------------------------------------------------Application: Bubble Point CalculationsGiven a liquid phase mixture, the bubble point is the temperature at which the first bubble of vapor appears, and the composition of that bubble. To estimate the bubble point temperature and composition, we invoke three assumptions: - Antoine's equation - provides an estimate of the pure component saturation pressure.
- Ideal gas - partial pressure of a component is equal to the vapor phase mole fraction times total pressure,
- Raoult's law - partial pressure of a component is equal to the liquid phase mole fraction times the saturation pressure,
The last two assumptions lead to a relationship which is the key to these vapor/liquid equilibrium problems. For bubble points, we know the pressure and liquid phase mole fractions . What's left to do is adjust the temperature until the vapor phase mole fractions add to one. Here we show how to formulate and solve for the bubble point of 50/50 mol% of benzene and toluene at 1 atm. % Use z to denote feed composition.P = 760;z_bnz = 0.5;z_tol = 0.5;% Using Raoult's law, solve for vapor phase mole fractions has a function% of temperaturey_bnz = @(T)z_bnz*bnz.Psat(T)/P;y_tol = @(T)z_tol*tol.Psat(T)/P;% The bubble point is the temperature at which the vapor phase mole% fractions sum to one.Tbub = fzero(@(T)y_bnz(T) + y_tol(T) - 1, 100);problem('Bubble Point Calculation')answer('z_bnz',z_bnz,'mole fraction');answer('z_tol',z_tol,'mole fraction');disp(sep);answer('Bubble Point Temperature',Tbub,'deg C');answer('Mole Fraction Benzene',y_bnz(Tbub),'mole fraction');answer('Mole Fraction Toluene',y_tol(Tbub),'mole fraction');% A graphical interpretation of the solution processT = 70:110;plot(T,y_bnz(T)+y_tol(T),T,y_bnz(T),T,y_tol(T));xlabel('Temperature [deg C]');ylabel('Mole Fraction');title('Benzene and Toluene vapor phase mole fractions');legend('y_{bnz}+y_{tol}','y_{bnz}','y_{tol}','Location','NW');ax = axis;hold on;plot([Tbub Tbub ax(1)],[ax(3) 1 1],'--');plot(Tbub,1,'.','Markersize',15);text(Tbub,1.05,sprintf(' T_{bubble} = %5.2f deg C',Tbub), ... 'HorizontalAlignment','center');plot(Tbub,y_bnz(Tbub),'g.','Markersize',15);plot([Tbub,ax(1)],y_bnz(Tbub)*[1 1],'g--');text(Tbub,y_bnz(Tbub),sprintf(' y_{bnz} = %5.3f',y_bnz(Tbub)));plot(Tbub,y_tol(Tbub),'r.','Markersize',15);plot([Tbub,ax(1)],y_tol(Tbub)*[1 1],'r--');text(Tbub,y_tol(Tbub),sprintf(' y_{tol} = %5.3f',y_tol(Tbub)));hold off;------------------------------------------------------------------Problem: Bubble Point Calculation------------------------------------------------------------------z_bnz = 0.5 [mole fraction]z_tol = 0.5 [mole fraction]------------------------------------------------------------------Bubble Point Temperature = 92.1 [deg C]Mole Fraction Benzene = 0.714 [mole fraction]Mole Fraction Toluene = 0.286 [mole fraction]Application: Dew Point CalculationGiven a vapor mixture at a fixed pressure, the dew point is the temperature at which the first drop of condensate appears. As for the bubble point computation, the key relationship is The we find the temperature for which the liquid phase mole fractions sum to one, % Use z to denote feed composition.P = 760;z_bnz = 0.5;z_tol = 0.5;% Solve for liquid phase mole fractions as a function of temperaturex_bnz = @(T)z_bnz*P./bnz.Psat(T);x_tol = @(T)z_tol*P./tol.Psat(T);% The bubble point is the temperature at which the vapor phase mole% fractions sum to one.Tdew = fzero(@(T)x_bnz(T) + x_tol(T) - 1, 100);problem('Dew Point Calculation')answer('z_bnz',z_bnz,'mole fraction');answer('z_tol',z_tol,'mole fraction');disp(sep);answer('Dew Point Temperature',Tdew,'deg C');answer('Mole Fraction Benzene',x_bnz(Tdew),'mole fraction');answer('Mole Fraction Benzene',x_tol(Tdew),'mole fraction');disp(sep);% A graphical interpretation of the solution processT = 70:110;plot(T,x_bnz(T)+x_tol(T),T,x_bnz(T),T,x_tol(T));xlabel('Temperature [deg C]');ylabel('Mole Fraction');title('Benzene and Toluene liquid phase mole fractions');legend('x_{bnz}+x_{tol}','x_{bnz}','x_{tol}','Location','NE');ax = axis;hold on;plot([Tdew Tdew ax(1)],[ax(3) 1 1],'--');plot(Tdew,1,'.','Markersize',15);text(Tdew,1.05,sprintf(' T_{dew} = %5.2f deg C',Tdew), ... 'HorizontalAlignment','center');plot(Tdew,x_bnz(Tdew),'g.','Markersize',15);plot([Tdew,ax(1)],x_bnz(Tdew)*[1 1],'g--');text(Tdew,x_bnz(Tdew),sprintf(' x_{bnz} = %5.3f',x_bnz(Tdew)));plot(Tdew,x_tol(Tdew),'r.','Markersize',15);plot([Tdew,ax(1)],x_tol(Tdew)*[1 1],'r--');text(Tdew,x_tol(Tdew),sprintf(' x_{tol} = %5.3f',x_tol(Tdew)));hold off;------------------------------------------------------------------Problem: Dew Point Calculation------------------------------------------------------------------z_bnz = 0.5 [mole fraction]z_tol = 0.5 [mole fraction]------------------------------------------------------------------Dew Point Temperature = 98.8 [deg C]Mole Fraction Benzene = 0.29 [mole fraction]Mole Fraction Benzene = 0.71 [mole fraction]------------------------------------------------------------------Application: Txy DiagramVapor/liquid equilibrium behavior of a binary mixture can be conveniently summarized in a Txy diagram. The Txy diagram depicts the solution to the equations The second of these equations comes from setting the partial pressures of the indivdual components equal to the total pressure. Solving the first equation for then substituting into the second equation % Fix total pressureP = 760;% Find the boiling points of the pure components. The Antoine class has a% method Tsat for doing this calculation.Tbnz = bnz.Tsat(P);Ttol = tol.Tsat(P);% A vector of temperatures at which to compute vapor and liquid phase% compositionsT = Tbnz:((Ttol-Tbnz)/100):Ttol;% Compute the liquid and vapor phase mole fractions for benzenex = (P-tol.Psat(T))./(bnz.Psat(T)-tol.Psat(T));y = x.*bnz.Psat(T)/P;% Plot the resultsplot(x,T,y,T);title('Txy diagram: benzene/toluene');xlabel('Mole fraction benzene');ylabel('Temperature [deg C]');legend('Dew Point','Bubble Point');% Show where the previous bubble and dew point calculations fit on the Txy% diagram.hold on;plot(x_bnz(Tdew),Tdew,'.','Markersize',15);plot(y_bnz(Tbub),Tbub,'g.','Markersize',15);plot([0.5 0.5],[Tbnz Ttol],'k--');plot([0.5 x_bnz(Tdew)],[Tdew Tdew],'b--');plot([0.5 y_bnz(Tbub)],[Tbub Tbub],'g--');hold off;Application: xy DiagramThe xy diagram provides much of the same information as the Txy diagram in format where the vapor phase composition is plotted as a function of the liquid phase composition. Temperature becomes a implicit parameter. % Reuse data for the Txy diagram to construct the xy diagramplot(x,y)xlabel('x_{bnz}');ylabel('y_{bnz}');axis('equal')axis('square');% Annotate xy diagram with temperaturehold on;for k = 1:10:length(T) plot(x(k),y(k),'.','Markersize',5); text(x(k),y(k),sprintf(' %4.1f',T(k)));endhold off% Show where previous bubble point and dew point calculations fit on this% plothold onplot([0.5 0.5 0],[0 y_bnz(Tbub) y_bnz(Tbub)],'g--');plot([0 x_bnz(Tdew) x_bnz(Tdew)],[0.5 0.5 0],'g--');hold offApplication: Isothermal FlashThe isothermal flash is a calculation coupling material balances with vapor/liquid equilibrium. Consider a flash unit with a molar feedrate and a vapor exit stream with flowrate and a liquid exit stream with molar flow . At steady state with component material balances where denotes the mole fraction of component i in the feed stream. From Raoult's law or where are the 'K-values'. Given the temperature, pressure, molar feedrate and feed compositions, the task is to compute the flowrates and composition of the exit streams. Traditionally, the first step in this calculation is to divide the first two equations by , and to define a vapor fraction phi. The second equation becomes Using Raoult's law to eliminate , solving each component balance for gives The liquid phase mole fractions must sum to one. Therefore, for a fixed temperature , we need to find such that Alternatively, one can compute the vapor phase mole fractions which must also sum to one Which of these equations to use to solve for ? While either can be made to work, it turns out that the difference between these equations has better numerical properties. The resulting equation is called the Rachford-Rice equation This is the equation we solve to find . % Feed and operating conditionsz_bnz = 0.5;z_tol = 0.5;% Operating ConditionsP = 760;T = 95;% K-valuesK_bnz = bnz.Psat(T)/P;K_tol = tol.Psat(T)/P;% Rachford-Rice equationr = @(phi) z_bnz*(K_bnz-1)./(1+phi*(K_bnz-1)) + ... z_tol*(K_tol-1)./(1+phi*(K_tol-1));% Solve for vapor phase fractionphi = fzero(r,[0,1]);% Show solutionproblem('Isothermal Flash')answer('Flash Temperature',T,'deg C');answer('Flash Pressure',P,'mmHg');answer('z_bnz',z_bnz,'mole fraction');answer('z_tol',z_tol,'mole fraction');disp(sep);answer('Vapor Phase Fraction',phi,'fraction');answer('y_bnz',K_bnz*z_bnz/(1+phi*(K_bnz-1)),'mole fraction');answer('y_tol',K_tol*z_tol/(1+phi*(K_tol-1)),'mole fraction');disp(sep);answer('Liquid Phase Fraction',1-phi,'fraction');answer('x_bnz',z_bnz/(1+phi*(K_bnz-1)),'mole fraction');answer('x_tol',z_tol/(1+phi*(K_tol-1)),'mole fraction');disp(sep);------------------------------------------------------------------Problem: Isothermal Flash------------------------------------------------------------------Flash Temperature = 95 [deg C]Flash Pressure = 760 [mmHg]z_bnz = 0.5 [mole fraction]z_tol = 0.5 [mole fraction]------------------------------------------------------------------Vapor Phase Fraction = 0.436 [fraction]y_bnz = 0.625 [mole fraction]y_tol = 0.375 [mole fraction]------------------------------------------------------------------Liquid Phase Fraction = 0.564 [fraction]x_bnz = 0.403 [mole fraction]x_tol = 0.597 [mole fraction]------------------------------------------------------------------
|
发表于 1970-1-1 08:00:00
显示全部楼层
IP卡
狗仔卡
发表于 1970-1-1 08:00:00
提升卡
置顶卡
喧嚣卡
变色卡
千斤顶
显身卡
成长值: 65875
成长值: 36850
成长值: 3940