View source

<source lang="a4c" tabwidth=4>
 REQUIRE "atoms.a4l";
 
 IMPORT "johnpye/extpy/extpy";
 
 IMPORT "vleplot";
 
 MODEL modifiedraoult;
 
 	(* first declare the variables you are working with *)
 
 	T IS_A temperature;
 
 	P, P1S, P2S IS_A pressure;
 
 	A IS_A factor;
 
 	gamma1, gamma2 IS_A positive_factor;
 
 	Pfac IS_A factor;
 	A1, B1, C1 IS_A positive_factor;
 	A2, B2, C2 IS_A positive_factor; x1, x2, y1, y2 IS_A fraction;
 
 	(* now write down the equations that describe the vapor liquid equilibrium *)
 	(* Antoine's equation for vapor pressure *) 
 
 	vp1: P1S/1000.0{Pa} = exp(A1 - B1/(T/1{K} - C1));
 	vp2: P2S/1000.0{Pa} = exp(A2 - B2/(T/1{K} - C2));
 
 	(* summation equations for mole fractions in vapor and liquid phases *)
 
 	sumx: x1 + x2 = 1.0;
 	sumy: y1 + y2 = 1.0;
 
 	(* the modified form of Raoult's Law - with a simple model for the activity coefficients *)
 
 	vle1: x1*gamma1*P1S = y1*P;
 	vle2: x2*gamma2*P2S = y2*P;
 
 	(* equations for the activity coefficients *)
 
 	activity1: gamma1 = exp(A*x2*x2);
 	activity2: gamma2 = exp(A*x1*x1);
 
 	eqA: A = 2.771 - 0.00523*T/1.0{K};
 
 	Pfac*100000.0 {Pa} = P;
 
 METHODS
 
     METHOD default_self;
         A1 := 16.59158;B1 := 3643.31;C1 := 33.424;
         A2 := 14.25326;B2 := 2665.54;C2 := 53.424;
     END default_self;
 
     METHOD specify;
     FIX A1, B1, C1, A2, B2, C2;
     FIX T, x1;
     END specify;
 
     METHOD values;
         T := 318.15 {K};
         x1 := 0.25;
     END values;
 
     METHOD on_load;
     RUN reset;
     RUN values;
     RUN default_self;
     END on_load;
 
         METHOD fancyplot;
     EXTERNAL vleplot(SELF);
         END fancyplot;
 
 END modifiedraoult; 
</source>

Here is the vleplot.py (most of the work by John)

<source lang="py" tabwidth=4>
 # python script for calculating Pressure versus y from T and x 
 # binary system, modified raoult's law, simple activity coefficient model
 
 import extpy
 from pylab import *
 from solverreporter import *
 
 def vleplot(self):
 
 	browser = extpy.getbrowser()
 	ioff()
 	figure()
 
 #
 # I have chosen three temperatures 
 #
 
 	for T in [320,340,370]:
 		self.T.setRealValue(T)
 #
 # collect the data for plotting in two sets of arrays (one for X, one for Y)
 # I have two sets here - one for P versus y and other for P versus x 
 #
 
 		XX1 = []
 		PP1 = []
 		XX2 = []
 		PP2 = []
 #
 # change x1 from 0 to 1.0 
 # there has to be a space between "in" and "["
 #
 
 		for x1 in [0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.99]:
 		    self.x1.setRealValue(x1)
 #
 # send the pair of values T x1 to the solver
 # and append the Pressure and y1 (from the solver) to the arrays 
 # the x's are also appended  
 		    try:
 				browser.sim.solve(browser.solver,SimpleSolverReporter(browser,message="T = %f, x1 = %f" % (T,x1)))
 				XX1.append(float(self.x1))
 				PP1.append(float(self.P))
 				XX2.append(float(self.y1))
 				PP2.append(float(self.P))
 		    except:
 		        browser.reporter.reportError('Failed to solve for x1 = %f' % x1)
 		        continue
 ## plot the data 
 		plot(XX1,PP1)
 		plot(XX2,PP2)
 		hold(1)
 	
 ##	legend()
 	ion()
 	show()
 
 extpy.registermethod(vleplot)
 #the above method can be called using "EXTERNAL vleplot(self)" in ASCEND.
 # if you want to see the azeotrope clearly, restrict the calculation to one
 # temperature
</source>

[[Category:Documentation]]