VLE examples

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;

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

# 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