Simple2Elegant

From ASCEND

This is an attempt at describing how models can be written for solutions in vapor liquid equilibrium.

We start by describing how to write MODELS and METHODS for example problem 10.1 from the book Chemical Engineering Thermodynamics by VanNess and Smith, 7th Edition, page 352. That is for a binary system of two components for which Raoult's Law is valid.

It is written in stages - the first stage being all variables defined as a "factor". Variables are later defined as "pressure" or "temperature" and so on. While ASCEND allows the modeler to write without much/any structure, as the models become more complex and as the need for their distribution becomes greater, paying attention to the art of model writing becomes important. Yet, when a new user starts with the ASCEND documentation, it can be difficult to read through the elegant models. Thus, this exercise.

This example ends with how Pxy or Txy plots can be generated python and the scripts written by John Pye.

A brief explanation of the variables. $P$ is the total system pressure, $x 1, x 2$ are the mole fractions of the two components in the liquid phase and $y 1, y 2$ are the mole fractions of the two components in the vapor phase and $P_1^S, P_2^S$ are the vapor pressures of components 1 and 2 (which we assume we can predict using Antoine's equation and $A_1, B_1, C_1, A_2, B_2~and~C_2$ are the Antoine's coefficients).

The relevant equations are (for a binary system)

$x_1 P_1^S = y_1 P$

$x_2 P_2^S = y_2 P$

$x 1 + x 2 = 1.0$

$y 1 + y 2 = 1.0$

$ln(P_1^S) = A_1 - \frac{B_1}{T + C_1}$

$ln(P_2^S) = A_2 - \frac{B_2}{T + C_2}$

Name of Calculation	What is Known	What we want to calculate
Bubble Point Pressure	T, x1	P, y1
Dew Point Pressure	T, y1	P, x1
Bubble Point Temperature	P, x1	T, y1
Dew Point Temperature	P, y1	T, x1

1 No Units, All variables declared as Factors
2 Using an exponential function, instead of natural log
3 Constants in a METHOD called values
4 Defining Units like Pressure and Temperature
5 Enter Temperature as Degrees C not in Kelvin
6 Define a MODEL for antoine so you can add more component values
7 Using the pygtk interface create Txy and Pxy plots
- 7.1 The python script that should be in ascdata or where ASCEND/pygtk can find it
8 Create your own models file/library file and include them when you want them
- 8.1 This file named mymodels.a4c has the model antoine and should be in ascdata (preferably) or share/ascend/models and you can add/edit
9 Additional Ideas

[edit] No Units, All variables declared as Factors

This is pg352version1.a4c

REQUIRE "atoms.a4l";

MODEL example101;
P1S, P2S, P, T, x1, x2, y1, y2 IS_A factor;

ln(P1S) = 14.2724 - 2945.47/(T + 224.0);
ln(P2S) = 14.2043 - 2972.64/(T + 209.0);
x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

METHODS
    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 75.0;
    END parta;
    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0;
        T := 300.0;
    END partb;
    METHOD doparta;
        RUN ClearAll;
    RUN parta;
    END doparta;
    METHOD dopartb;
        RUN ClearAll;
    RUN partb;
    END dopartb;

END example101;

[edit] Using an exponential function, instead of natural log

This is pg352version2.a4c

REQUIRE "atoms.a4l";

MODEL example101;

P1S, P2S, P, T, x1, x2, y1, y2 IS_A factor;

P1S = exp(14.2724 - 2945.47/(T + 224.0));
P2S = exp(14.2043 - 2972.64/(T + 209.0));

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

METHODS
    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 75.0;
    END parta;
    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0;
        T := 300.0;
    END partb;
    METHOD doparta;
        RUN ClearAll;
    RUN parta;
    END doparta;
    METHOD dopartb;
        RUN ClearAll;
    RUN partb;
    END dopartb;

END example101;

[edit] Constants in a METHOD called values

Instead of hardwiring the antoine constants, we put them in a METHOD called values This is pg352version3.a4c

REQUIRE "atoms.a4l";

MODEL example101;

A1, B1, C1, A2, B2, C2 IS_A factor;
P1S, P2S, P, T, x1, x2, y1, y2 IS_A factor;

ln(P1S) = A1 - B1/(T + C1);
ln(P2S) = A2 - B2/(T + C2);

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

METHODS
    METHOD values;
    A1 := 14.2724; B1 := 2945.47; C1 := 224.0;
    A2 := 14.2043; B2 := 2972.64; C2 := 209.0;
    FIX A1; FIX B1; FIX C1; FIX A2; FIX B2; FIX C2;
    END values;
    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 75.0;
    END parta;
    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0;
        T := 300.0;
    END partb;
    METHOD doparta;
        RUN ClearAll;
    RUN values;
    RUN parta;
    END doparta;
    METHOD dopartb;
        RUN ClearAll;
    RUN values;
    RUN partb;
    END dopartb;

END example101;

[edit] Defining Units like Pressure and Temperature

Note how we had to transform the vapor pressure equations. ASCEND use SI units internally This is pg352version4.a4c

REQUIRE "atoms.a4l";

MODEL example101;

A1, B1, C1, A2, B2, C2 IS_A factor;
P1S, P2S, P IS_A pressure;
T IS_A temperature;
x1, x2, y1, y2 IS_A fraction;

P1S/100000.0{Pa} = exp(A1 - B1/(T/1.0{K} - 273.15 + C1));
P2S/100000.0{Pa} = exp(A2 - B2/(T/1.0{K} - 273.15 + C2));

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

METHODS
    METHOD default_self;
    END default_self;
    METHOD values;
    A1 := 14.2724; B1 := 2945.47; C1 := 224.0;
    A2 := 14.2043; B2 := 2972.64; C2 := 209.0;
    FIX A1; FIX B1; FIX C1; FIX A2; FIX B2; FIX C2;
    END values;

    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 348.15 {K};
    END parta;
    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0 {kPa};
    END partb;
    METHOD doparta;
        RUN ClearAll;
    RUN values;
    RUN parta;
    END doparta;
    METHOD dopartb;
        RUN ClearAll;
    RUN values;
    RUN partb;
    END dopartb;

END example101;

[edit] Enter Temperature as Degrees C not in Kelvin

This is pg352version5.a4c

REQUIRE "atoms.a4l";

MODEL example101;

A1, B1, C1, A2, B2, C2 IS_A factor;
P1S, P2S, P IS_A pressure;
T IS_A temperature;
x1, x2, y1, y2 IS_A fraction;
T_degC IS_A factor;

P1S/1000.0{Pa} = exp(A1 - B1/(T/1.0{K} - 273.15 + C1));
P2S/1000.0{Pa} = exp(A2 - B2/(T/1.0{K} - 273.15 + C2));

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

T_degC = T/1{K} - 273.15;

METHODS
    METHOD default_self;
    END default_self;
    METHOD values;
    A1 := 14.2724; B1 := 2945.47; C1 := 224.0;
    A2 := 14.2043; B2 := 2972.64; C2 := 209.0;
    FIX A1; FIX B1; FIX C1; FIX A2; FIX B2; FIX C2;
    END values;

    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 348.15 {K};
    END parta;
    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0 {kPa};
    END partb;
    METHOD doparta;
        RUN ClearAll;
    RUN values;
    RUN parta;
    END doparta;
    METHOD dopartb;
        RUN ClearAll;
    RUN values;
    RUN partb;
    END dopartb;

END example101;

[edit] Define a MODEL for antoine so you can add more component values

This is pg352version6.a4c

REQUIRE "atoms.a4l";

MODEL antoine(nc WILL_BE symbol_constant;);

A, B, C IS_A factor;

SELECT(nc)

CASE 'acetonitrile':A = 14.2724; B = 2945.47; C = 224.070;
CASE 'nitromethane':A = 14.2043; B = 2972.64; C = 209.000;

END SELECT;

END antoine;

MODEL example101;

A1, B1, C1, A2, B2, C2 IS_A factor;
P1S, P2S, P IS_A pressure;
T IS_A temperature;
x1, x2, y1, y2 IS_A fraction;
T_degC IS_A factor;

nc1, nc2 IS_A symbol_constant;

nc1 :== 'acetonitrile';
nc2 :== 'nitromethane';

mync1 IS_A antoine(nc1);
mync2 IS_A antoine(nc2);

A1 = mync1.A;
B1 = mync1.B;
C1 = mync1.C;
A2 = mync2.A;
B2 = mync2.B;
C2 = mync2.C;

P1S/1000.0{Pa} = exp(A1 - B1/(T/1.0{K} - 273.15 + C1));
P2S/1000.0{Pa} = exp(A2 - B2/(T/1.0{K} - 273.15 + C2));

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

T_degC = T/1{K} - 273.15;

METHODS

    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 348.15 {K};
    END parta;
    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0 {kPa};
    END partb;
    METHOD doparta;
        RUN ClearAll;
    RUN parta;
    END doparta;
    METHOD dopartb;
        RUN ClearAll;
    RUN partb;
    END dopartb;

END example101;

[edit] Using the pygtk interface create Txy and Pxy plots

This is pg352version7.a4c

REQUIRE "atoms.a4l";
IMPORT "johnpye/extpy/extpy";
IMPORT "vleplots";

MODEL antoine(nc WILL_BE symbol_constant;);

A, B, C IS_A factor;

SELECT(nc)

CASE 'acetonitrile':A = 14.2724; B = 2945.47; C = 224.070;
CASE 'nitromethane':A = 14.2043; B = 2972.64; C = 209.000;

END SELECT;

END antoine;

MODEL example101;

A1, B1, C1, A2, B2, C2 IS_A factor;
P1S, P2S, P IS_A pressure;
T IS_A temperature;
x1, x2, y1, y2 IS_A fraction;
T_degC IS_A factor;
TdegC IS_A factor;

nc1, nc2 IS_A symbol_constant;

nc1 :== 'acetonitrile';
nc2 :== 'nitromethane';

mync1 IS_A antoine(nc1);
mync2 IS_A antoine(nc2);

A1 = mync1.A;
B1 = mync1.B;
C1 = mync1.C;
A2 = mync2.A;
B2 = mync2.B;
C2 = mync2.C;

P1S/1000.0{Pa} = exp(A1 - B1/(T/1.0{K} - 273.15 + C1));
P2S/1000.0{Pa} = exp(A2 - B2/(T/1.0{K} - 273.15 + C2));

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

T_degC = T/1{K} - 273.15;
TdegC = T_degC;

METHODS

    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 348.15 {K};
    END parta;

    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0 {kPa};
    END partb;

    METHOD doparta;
        RUN ClearAll;
    RUN parta;
    END doparta;

    METHOD dopartb;
        RUN ClearAll;
    RUN partb;
    END dopartb;

        METHOD generatepxyplot;
    RUN ClearAll;
    RUN parta;
    RUN doparta;
        EXTERNAL pxyplot(SELF);
        END generatepxyplot;

        METHOD generatetxyplot;
    RUN ClearAll;
    RUN partb;
    RUN dopartb;
        EXTERNAL txyplot(SELF);
        END generatetxyplot;

END example101;

[edit] The python script that should be in ascdata or where ASCEND/pygtk can find it

Script written by John Pye, minor changes by Krishnan Chittur. There are two methods txy and pxyplot This is vleplots.py

import extpy
from pylab import *
from solverreporter import *
def txyplot(self):
    browser = extpy.getbrowser()
    ioff()
    figure()
    for P in [70000.0]:
        self.P.setRealValue(P)
        XX1 = []
        TT1 = []
        XX2 = []
        TT2 = []
        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)
            try:
                browser.sim.solve(browser.solver,SimpleSolverReporter(browser,message="P = %f, x1 = %f" % (P,x1)))
                XX1.append(float(self.x1))
                TT1.append(float(self.TdegC))
                XX2.append(float(self.y1))
                TT2.append(float(self.TdegC))
            except:
                browser.reporter.reportError('Failed to solve for x1 = %f' % x1)
                continue
        plot(XX1,TT1)
        plot(XX2,TT2)
        hold(1)
    ion()
    show()

extpy.registermethod(txyplot)

def pxyplot(self):
    browser = extpy.getbrowser()
    ioff()
    figure()
    for T in [340]:
        self.T.setRealValue(T)
        XX1 = []
        PP1 = []
        XX2 = []
        PP2 = []
        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)
            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(XX1,PP1)
        plot(XX2,PP2)
        hold(1)
    ion()
    show()

extpy.registermethod(pxyplot)

[edit] Create your own models file/library file and include them when you want them

This is pg352version8.a4c

REQUIRE "atoms.a4l";
REQUIRE "mymodels.a4c";
IMPORT "johnpye/extpy/extpy";
IMPORT "vleplots";

MODEL example101;

A1, B1, C1, A2, B2, C2 IS_A factor;
P1S, P2S, P IS_A pressure;
T IS_A temperature;
x1, x2, y1, y2 IS_A fraction;
T_degC IS_A factor;
TdegC IS_A factor;

nc1, nc2 IS_A symbol_constant;

nc1 :== 'acetonitrile';
nc2 :== 'nitromethane';

mync1 IS_A antoine(nc1);
mync2 IS_A antoine(nc2);

A1 = mync1.A;
B1 = mync1.B;
C1 = mync1.C;
A2 = mync2.A;
B2 = mync2.B;
C2 = mync2.C;

P1S/1000.0{Pa} = exp(A1 - B1/(T/1.0{K} - 273.15 + C1));
P2S/1000.0{Pa} = exp(A2 - B2/(T/1.0{K} - 273.15 + C2));

x1 + x2 = 1.0;
y1 + y2 = 1.0;

x1*P1S = y1*P;
x2*P2S = y2*P;

T_degC = T/1{K} - 273.15;
TdegC = T_degC;

METHODS

    METHOD parta;
    FIX T; FIX x1;
    x1 := 0.6;
    T := 348.15 {K};
    END parta;

    METHOD partb;
    FIX P; FIX x1;
    x1 := 0.1;
    P := 70.0 {kPa};
    END partb;

    METHOD doparta;
        RUN ClearAll;
    RUN parta;
    END doparta;

    METHOD dopartb;
        RUN ClearAll;
    RUN partb;
    END dopartb;

        METHOD generatepxyplot;
    RUN ClearAll;
    RUN parta;
    RUN doparta;
        EXTERNAL pxyplot(SELF);
        END generatepxyplot;

        METHOD generatetxyplot;
    RUN ClearAll;
    RUN partb;
    RUN dopartb;
        EXTERNAL txyplot(SELF);
        END generatetxyplot;

END example101;

[edit] This file named mymodels.a4c has the model antoine and should be in ascdata (preferably) or share/ascend/models and you can add/edit

This is mymodels.a4c

REQUIRE "atoms.a4l";

MODEL antoine(nc WILL_BE symbol_constant;);

A, B, C IS_A factor;

SELECT(nc)

CASE 'acetonitrile':A = 14.2724; B = 2945.47; C = 224.070;
CASE 'nitromethane':A = 14.2043; B = 2972.64; C = 209.000;

END SELECT;

END antoine;

[edit] Additional Ideas

Extend this to multiple components (more than 2) and create a components file for antoine coefficients

Retrieved from "http://ascendwiki.cheme.cmu.edu/Simple2Elegant"

Category: Documentation

Simple2Elegant

From ASCEND

Contents

[edit] No Units, All variables declared as Factors

[edit] Using an exponential function, instead of natural log

[edit] Constants in a METHOD called values

[edit] Defining Units like Pressure and Temperature

[edit] Enter Temperature as Degrees C not in Kelvin

[edit] Define a MODEL for antoine so you can add more component values

[edit] Using the pygtk interface create Txy and Pxy plots

[edit] The python script that should be in ascdata or where ASCEND/pygtk can find it

[edit] Create your own models file/library file and include them when you want them

[edit] This file named mymodels.a4c has the model antoine and should be in ascdata (preferably) or share/ascend/models and you can add/edit

[edit] Additional Ideas

Views

Personal tools

Navigation

tools

wiki stuff

Search

Toolbox