978-1107135116 Chapter 3 Part 4

Nonlinear equations 179

53 xlabel(‘$x$’,‘FontSize’,14)

54 ylabel(‘$y$’,‘FontSize’,14)

55 title(‘Solution to s13c6p3’,‘FontSize’,14)

56 saveas(g,‘s13c6p123 solution figure3.eps‘)

67 J(1,2) = -y;

68 J(2,1) = P;

69 J(2,2) = -(1-y);

70 out = J;

The two output files are:

0 0.2 0.4 0.6 0.8 1

800

900

1000

1100

1200

1300

1400

1500

x, y

P(mm H g)

Solution t o s 13c 6p3

0 0.2 0.4 0.6 0.8 1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

y

Solution to s13c 6p3

180 Nonlinear equations

Note on the thermodynamics problem. One nice thing about the one param-

eter model is that you can solve for the azeotrope without making a nonlinear

problem. The root of f(x1) is

2

problem 3.

1function out = s13c6p1 fast

2clc

3close all

4set(0,‘defaulttextinterpreter’,‘latex’)

15 xlabel(‘$x {\mathrm{azeo}}$’,‘FontSize’,14)

16 ylabel(‘$P {\mathrm{azeo}}$’,‘FontSize’,14)

17 title(‘Solution to s13c6p1’,‘FontSize’,14)

18 saveas(h,‘s13c6p1 fast plot.eps’,‘psc2’)

1function out = s13c6p3 fast

11 npts = 200;

12 x = linspace(0,1,npts); %values of x1

13 P = x.*exp(A*(1-x).ˆ2)*P1sat + (1-x).*exp(A*x.ˆ2)*P2sat;

14 y = x.*exp(A*(1-x).ˆ2)*P1sat./P;

15

26 ylabel(‘$y$’,‘FontSize’,14)

27 title(‘Solution to s13c6p3’,‘FontSize’,14)

28 saveas(g,‘s13c6p3 fast plot2.eps’,‘psc2’)

Note that the solution here only works for the one-parameter model. While this

is the best way to solve the problem for a thermodynamics class, it gives you no

(3.31) (a) The files for this problem are contained in the folder solutions/s15c6p1 matlab.

At the azeotrope, x=y. If we divide the two equilibrium expressions we get

γ1

=Psat

2

1−β2x3

1−ln Psat

Psat

1

The derivative we need is

f′=−2α1(1 −x1)−3β1(1 −x1)2−2α2x1−3β2x2

8% beta1 = 0.1;

9% beta2 = 0.05;

10 % P1sat = 1100;

11 % P2sat = 800;

12

23 %put all the parameters into a vector for easy passing

24 p(1,1) = alpha1;

25 p(2,1) = alpha2;

26 p(3,1) = beta1;

27 p(4,1) = beta2;

38 f = getf(x,p);

39 k=k+1;

40 if k>20

41 fprintf(‘Did not converge.\n’)

42 end

53 y=x*P1sat*gamma1/P;

54 fprintf(‘y = %8.6f\n’,y)

55

56 %Check that this is the maximum pressure

Nonlinear equations 183

67 beta1 = p(3);

68 beta2 = p(4);

69 P1sat = p(5);

70 P2sat = p(6);

71 out = alpha1*(1-x)ˆ2 + beta1*(1-x)ˆ3 – alpha2*xˆ2 – …

3*beta2*xˆ2;

The code is a bit more complicated than necessary, but I have set it up so that

it will be easier to use for later problems. The result is xazeo =0.744149 at

Pazeo =1138.5037 mm Hg.

(b) The files to solve this problem are contained in the folder solutions/s15c6p2 matlab.

x=





γ1

γ2

x1

y1







0 (1 −x1)Psat

2−γ2Psat

2P

The Matlab file to implement Newton-Raphson using the azeotrope as the initial

guess is:

184 Nonlinear equations

1function s15h6p2

2clc

13 P2sat = 800;

14 P = 1000; %junk value for now to put in place of pressure

15

16 %put all the parameters into a vector for easy passing

17 p(1,1) = alpha1;

24 %get the azeotope values from Problem #1

25 [x,P] = get azeotrope(p);

26 p(7,1) = P;

27

28 %set tolerance for Newton-Raphson

39

40 function [x,flag] = nonlinear newton raphson(x,tol,p)

41 % x = initial guess for x

42 % tol = criteria for convergence

43 % requires function R = getR(x)

54 J = getJ(x,p); %compute Jacobian

55 del = -J\R; %Newton Raphson

56 x = x + del; %update x

Nonlinear equations 185

57 k=k+1;%update counter

68 end

69

70

71 function out = getR(x,p)

72 %unpack the parameters for clarity

83 gamma2 = x(2);

84 x1 = x(3);

85 y1 = x(4);

86

87 %write equations

98 alpha1 = p(1);

99 alpha2 = p(2);

100 beta1 = p(3);

101 beta2 = p(4);

102 P1sat = p(5);

186 Nonlinear equations

113 J = zeros(4);

114 J(1,1) = 1/gamma1;

115 J(1,3) = 2*alpha1*(1-x1) + 3*beta1*(1-x1)ˆ2;

116 J(2,2) = 1/gamma2;

127 function [xvec,P] = get azeotrope(p)

128 %unpack parameters

129 alpha1 = p(1);

130 alpha2 = p(2);

131 beta1 = p(3);

142 x = x – f/df;

143 f = getf(x,p);

144 k=k+1;

145 if k>20

146 fprintf(‘Did not converge.\n’)

157 %pack into Newton-Raphson vector

158 xvec(1,1) = gamma1;

159 xvec(2,1) = gamma2;

160 xvec(3,1) = x;

161 xvec(4,1) = x;

beta2*xˆ3 – log(P2sat/P1sat);

171

172 function out = getdf(x,p)

173 alpha1 = p(1);

174 alpha2 = p(2);

even if you gave a good guess, you will get a bad solution for δdue to a high

condition number, which means the new value of xmay not be better than the

previous guess. You can see how Newton-Raphson blows up when you set the

tolerance lower so that the while loop executes. There are tricks to handling

poorly conditioned matrices, but it seems unnecessary here since the residual is

2) or towards x1=1. Another option is to start from

cases very close to the pure species, where we know that the activity coeﬃcients

are γ≈1 and the pressure is almost the saturation pressure. If you use this second

approach, you need to increase the pressure until you get near the azeotrope. I did

both approaches in this file to show you how it is done. In either case, be careful

10 alpha2 = A – 3*B;

11 beta1 = -4*B;

12 beta2 = 4*B;

188 Nonlinear equations

13 P1sat = 1100;

24 %get the azeotope values from Problem #1

25 [xazeo,Pazeo] = get azeotrope(p);

26 P = Pazeo;

27 p(7,1) = P;

28

38 x1 = xazeo-0.02;

39 y1 = xazeo-0.05;

40 gamma1 = exp(alpha1*(1-x1)ˆ2 + beta1*(1-x1)ˆ3);

41 gamma2 = exp(alpha2*x1ˆ2 + beta2*x1ˆ3);

42 x(1,1) = gamma1;

53 xleft(i) = x(3);

54 yleft(i) = x(4);

55 Pleft(i) = P;

56

57 fprintf(‘The solutions for P = %8.2f is\n’,P)

Nonlinear equations 189

68 fprintf(‘***************************************\n’)

69

70

71 %now start from almost pure species 2 and go up

82 npts = floor((Pazeo-Pmin)/dP);

83 xright = zeros(npts,1);

84 yright = zeros(npts,1);

85 Pright = zeros(npts,1);

86

97 fprintf(‘The solutions for P = %8.2f is\n’,P)

98 fprintf(‘ gamma1 = %8.6f \n’,x(1))

99 fprintf(‘ gamma2 = %8.6f \n’,x(2))

100 fprintf(‘ x1 = %8.6f \n’,x(3))

101 fprintf(‘ y1 = %8.6f \n’,x(4))

111 ylabel(‘Pressure (mm Hg)’,‘FontSize’,14)

112 saveas(h,‘s15h6p3 solution figure.eps’,‘psc2’)

113

114

115 function [x,flag] = nonlinear newton raphson(x,tol,p)

125 while norm(R) >tol

126 J = getJ(x,p); %compute Jacobian

127 del = -J\R; %Newton Raphson

128 x = x + del; %update x

129 k=k+1;%update counter

140 end

141

142 function out = getR(x,p)

143 %unpack the parameters for clarity

144 alpha1 = p(1);

155 x1 = x(3);

156 y1 = x(4);

157

158 %write equations

159 R = zeros(4,1);

170 alpha2 = p(2);

171 beta1 = p(3);

172 beta2 = p(4);

173 P1sat = p(5);

174 P2sat = p(6);

184 J = zeros(4);

185 J(1,1) = 1/gamma1;

186 J(1,3) = 2*alpha1*(1-x1) + 3*beta1*(1-x1)ˆ2;

187 J(2,2) = 1/gamma2;

188 J(2,3) = -2*alpha2*x1 – 3*beta2*x1ˆ2;

199 %unpack parameters

200 alpha1 = p(1);

201 alpha2 = p(2);

202 beta1 = p(3);

203 beta2 = p(4);

214 f = getf(x,p);

215 k=k+1;

216 if k>20

217 fprintf(‘Did not converge.\n’)

218 end

229 alpha1 = p(1);

230 alpha2 = p(2);

231 beta1 = p(3);

232 beta2 = p(4);

233 P1sat = p(5);

242 P1sat = p(5);

243 P2sat = p(6);

244 out = -2*alpha1*(1-x) – 3*beta1*(1-x)ˆ2 – 2*alpha2*x – …

3*beta2*xˆ2;

The output file is:

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

800

850

900

950

1000

1050

1100

1150

1200

1250

1300

C on c e nt rat i on x1( bl ue ) or y1( re d )

P r e ss ur e ( mm H g)

I made the azeotrope as a green xso that it is easy to see.

(3.32) The files for this problem are contained in the folder s12c5p12345 matlab.

The Matlab file that solves this entire problem is:

1function s12c5p12345

2clc

3close all

13

14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

15 %%% QUESTION #1

16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

17 fprintf(‘\n\n’)

27 V guess = 0.5;

28 V = find V(P,V guess,T,a,b); %use Newton’s method to find V

29 fprintf(‘V = %8.6f at P = %3.2f MPa and T = %6.2f K\n’,V,P,T)

30

31

42 h = plot vdw(T,a,b,0.21,3,0,6);

43 saveas(h,‘s12c5p12345 solution figure2a.eps’,‘psc2’)

44

45 %use a decent guess

46 V = 0.23;

57 fprintf(‘PROBLEM #3 \n’)

58 fprintf(‘*****************************************************\n’)

59

60 %use a bad guess

61 V = 0.6;

72 fprintf(‘\n\n’)

73 fprintf(‘*****************************************************\n’)

74 fprintf(‘PROBLEM #4 \n’)

194 Nonlinear equations

75 fprintf(‘*****************************************************\n’)

86 %compute the values of Vl, Vg and Psat at this T

87 [Vl,Vg,Psat] = equilibrium(P,T,Vl,Vg,a,b);

88 %plot the isotherm

89 h = plot isotherm(T,Vl,Vg,Psat,a,b);

90 saveas(h,‘s12c5p12345 solution figure4b.eps’,‘psc2’)

101 %use the values from previous steps to make the phase envelope

102 fprintf(‘Generating phase envelope. Will output results at …

each step.\n’)

103 [g,h] = phase envelope(Pc,Vc,Tc,T,Psat,Vl,Vg,a,b);

104 saveas(h,‘s12c5p12345 solution figure5a.eps’,‘psc2’)

114 Tc = (8*a)/(27*b*R);

115 Pc = (a/bˆ2)/27;

116

117 %display the results to the screen

118 fprintf(‘For a = %6.5f MPa Lˆ2/molˆ2 and b = %5.4f L/mol, \n’,a,b)

Nonlinear equations 195

128 for i = 1:npts

129 P(i) = vdw(V(i),T,a,b);

130 end

131 h = figure;

the molar

142 %volume in L/mol and the temperature in K, a in in units of …

MPa*Lˆ2/molˆ2

143 %and b in units of L/mol

144 R = 0.008314; %ideal gas constant in MPa*L/mol K

153 tol = 1e-8;

154 fVeval = fvdw(P,V,T,a,b); %value of function f(V)

155 k = 0; %counter

156 while abs(fVeval) >tol

157 V = V – fVeval/dfvdw(V,T,a,b);

167 out = V;

168

169 function [h,kcount,errcount] = find V plot(P,V,T,a,b)

170 %executes Newton’s method to find a value of V in L/mol when …

the vdw

the figure too.

176 Vinit = V;

177 tol = 1e-8;

178 fVeval = fvdw(P,V,T,a,b); %value of function f(V)

179 k = 0; %counter

190 fVeval = fvdw(P,V,T,a,b);

191 k = k+1;

192 kcount(k+1) = k;

193 Vcount(k+1) = V;

194 errcount(k+1) = abs(fVeval);

204 h=figure;

205 plot(kcount,Vcount,‘-ob’)

206 xlabel(‘Newton method iteration, $k$’,‘FontSize’,14)

207 ylabel(‘Value of molar volume, $\underline{V}ˆ{k}$…

(L/mol)’,‘FontSize’,14)

216 kscale(npts) = kcount(npts);

217 escale(npts) = errcount(npts);

218 for i = npts-1:-1:1

219 kscale(i) = kcount(i);

220 escale(i) = errcount(i+1)ˆ(1/2);

229 title(top,‘FontSize’,14)

230

231

232 function out = fvdw(P,V,T,a,b)

233 %computes the value of the vdw equation in the form f(V) = 0 …

242 out = -R*T/(V-b)ˆ2+2*a/Vˆ3;

243

244 function out = plot Psat guess(P,T,a,b)

245 %make a plot of the initial guess. returns figure handle for …

saving

256 axis([Vmin,Vmax,0,3])

257 xlabel(‘$\underline{V}$ (L/mol)’,‘FontSize’,14)

258 ylabel(‘$P$ (MPa)’,‘FontSize’,14)

259 top = strcat(‘Guess of Psat for $T = $’,num2str(T),‘K’);

260 title(top,‘FontSize’,14)

270 Vg = find V(Psat,Vg,T,a,b);

271 fGeval = fGdep(Psat,Vl,Vg,T,a,b);

272

273 %start iterating

274 k = 1;

282 %update the volumes

283 Vl = find V(Psat,Vl,T,a,b);

284 Vg = find V(Psat,Vg,T,a,b);

285 %compute the new dG

286 fGeval = fGdep(Psat,Vl,Vg,T,a,b);

296

297 function out = fGdep(P,Vl,Vg,T,a,b)

298 %computes the difference in the departure Gibbs free energy …

for the vdw

299 %equation with P in MPa, V in L/mol, T in K and returns dG in …

309 out = (Vl-Vg);

310

311 function out = plot isotherm(T,Vl,Vg,Psat,a,b)

312 %create an isotherm below the critical point. returns handle …

for plotting