# 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;

RUN reset;
RUN values;

RUN default_self;

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 &quot;in&quot; and &quot;[&quot;
#

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=&quot;T =&nbsp;%f, x1 =&nbsp;%f&quot; % (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 =&nbsp;%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 &quot;EXTERNAL vleplot(self)&quot; in ASCEND.
# if you want to see the azeotrope clearly, restrict the calculation to one

# temperature
```