Figure 3.6:

Gibbs energy contours for the pentane reactions as a function of the two reaction extents.

Code for Figure 3.6

Text of the GNU GPL.

main.m


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

% Revised 7/24/2018

p = struct();
p.deltag1 = -3.72; % kcal/mol
p.deltag2 = -4.49; % kcal/mol
p.T       = 400;   % K
p.P       = 2.5;   % atm
p.R       = 1.987; % cal/mol/K
p.K1 = exp(-p.deltag1*1000/(p.R*p.T));
p.K2 = exp(-p.deltag2*1000/(p.R*p.T));




% first optimize free energy
% min_z G(z) s.t. x_1,x_2 >0,  x_1 + x_2 < 0.5

p.Gval = 0;
p.w0 = [0.2; 0.2];
p.A = [1, 1];
p.b = 0.5;
[w, obj, info] = fmincon(@(w) G(w,p), p.w0, p.A, p.b, [], [], [0;0], [0.5;0.5]);
w

if (info ~= 1)
    info, w, obj
end

% next prepare contour plot

Glevels = [0, -1, -2, -2.5, -2.5, -2.53, -2.55 -2.559];
nlev = length(Glevels);
alpha = [0.05, 0.05,   0.075, 0.05, 0.05, 0.05, 0.04, 0.02];
x0vec = [1e-3, 1e-3, 1e-3, .133, .133, .133, .133, .133; ...
          .02,  .15,  .35,   .3,  .3,   .3,    .3,   .3];
delx  = [0.01, 0.01, 0.01, 0.01, -0.01, 0, 0, 0; ...
            0     0,    0,    0,     0, 0, 0, 0];
sfinal = [9, 9, 9, 9, 9, 9, 9, 5];

ztable(1:nlev) = {NaN(2,1)};
for i = 1:nlev
    p.Gval = Glevels(i);

    %  find a point on the contour corresponding to s=0
    ext0 = x0vec(:,i);
    unk = [2];
    x0 = ext0(unk);
    [x, fval, info] = fsolve(@(x) Gcontour(x,p,ext0,unk), x0);
    fsolve_failed = info <= 0;
    infovec(i) = info;
    fcnvec(:,i) = Gcontour(x,p,ext0,unk);
    ext = ext0;
    ext(unk) = x;
    zlast = ext;
    extinit(:,i) = ext;

    % continuation in arclength
    dels = alpha(i);
    svec = dels:dels:sfinal(i);
    narc =length(svec);
    z0 = ext + delx(:,i);
    unk = [1 2];

    for j = 1:narc
        s = svec(j);
        [z, fval, info] = fsolve(@(z) inflate(z,zlast,dels,p,ext0,unk),z0);
        fsolve_ok = info > 0;
        % stop when contour completed

        if (z(1) > 0 && z(2) > 0 && 0.5-(z(1)+z(2)) > 0 && fsolve_ok)
            ztable{i}(:,j) = z;
            % project initial guess along the arc for next point on contour
            z0 = 2*z - zlast;
            zlast = z;
        else
            break
        end
    end

          ztable(i) = {[ext, ztable{i}]'};
    end

feasbound = [0, 0.5; 0.5, 0];  % store boundary of feasible region
save contour_pentane.dat ztable feasbound;


%plot results
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
    hold on;
    plot([0.5,0], [0, 0.5]); % one leg of feasible region

    for i = 1:nlev
        plot(ztable{i}(:,1), ztable{i}(:,2));
    end

    % TITLE
end % PLOTTING

Gcontour.m


1
2
3
4
function retval = Gcontour(var,p,ext0,unk)
    w = ext0;
    w(unk) = var;
    retval = G(w,p);

G.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
function retval = G(w,p)
    x = w(1);
    y = w(2);
    z = 0.5 - x - y;
    pp = 0.5 + z;

    if (x == 0)
        xlogx = 0;
    else
        xlogx = abs(x)*log(abs(x/pp));
    end

    if (y == 0)
        ylogy = 0;
    else
        ylogy = abs(y)*log(abs(y/pp));
    end

    if (z == 0)
        zlogz = 0;
    else
        zlogz = abs(z)*log(abs(z/pp));
    end

    retval = -p.Gval -(x*log(p.K1) + y*log(p.K2)) + pp*log(p.P) + ...
    2*zlogz + xlogx + ylogy;

inflate.m


1
2
function retval = inflate(z,zlast,dels,p,ext0,unk)
    retval = [ Gcontour(z,p,ext0,unk) ; norm(z-zlast) - dels^2 ];