Cutlip Cutli p and Shacham: Probl Problem em Solving in Chemi Chemical cal and Biochem Biochemical ical Engineering Engineering
Chapter 5
Problem Solving with MATLAB
Cheng-Liang Chen
PSE
LABORATORY
Department of Chemical Engineering National TAIWAN University
Chen CL
1
Molar Volume and Compressibility from Redlich-Kwong Equation Analytical ytical solution solution of the cubic Redlich-Kwo Redlich-Kwong ng equation for Concepts Utilized: Anal compressibility factor and calculation of the molar volume at various reduced temperature and pressure values. Numerical Methods: Solution of a set of explicit equations. Problem Statement:
The Redlich-Kwong equation of state is given by
P =
RT (V b)
a √ − − V V ((V + b) T
5/2
R2T c a = 0.42747 P c RT c b = 0.08664 P c
(1)
Chen CL
1
Molar Volume and Compressibility from Redlich-Kwong Equation Analytical ytical solution solution of the cubic Redlich-Kwo Redlich-Kwong ng equation for Concepts Utilized: Anal compressibility factor and calculation of the molar volume at various reduced temperature and pressure values. Numerical Methods: Solution of a set of explicit equations. Problem Statement:
The Redlich-Kwong equation of state is given by
P =
RT (V b)
a √ − − V V ((V + b) T
5/2
R2T c a = 0.42747 P c RT c b = 0.08664 P c
(1)
Chen CL
2
where P = pre pressure ssure in atm V = mola molarr volume volume in liters/g-m liters/g-mol ol T = temperat temperature ure in K R = gas const constant ant (R = 0.08206 08206 (atm-liter/g-mol-K)) (atm-liter/g-mol-K)) T e = critica criticall temperature temperature in K P e = critica criticall pressure pressure in atm atm The compressibility factor is given by z=
P V RT
Equation (1) can be written, after considerable algebra, in terms of the compressibility factor as a cubic equation (see Seader and Henley) f ((z ) = z 3 f
2
− z − qz − r = 0
where r = A2B, 2
A
= 0.42747
(5)
q = B 2 + B
P r 5/2 T r
B
= 0.0866
2
−A
P r T r
in which P r is the reduced pressure (P ( P /P c) and T r is the reduced temperature (T /T c). Equation (5) can be solved analytically for three roots. Some of these
Chen CL
3
roots are complex. Considering only the real roots, the sequence of calculations involves the steps 3
C = f =
f g 2 + 3 2 3q 1 g = 3
− −
−27r − 9q − 2 27
If C > 0, there is one real solution for z given by z = D + E + 1/3 D =
√ (−g/2 + C ) /
1 3
E =
√ (−g/2 − C ) /
1 3
If C < 0, there are three real solutions
−
f φ 2π(k 1) 1 zk = 2 cos + + , 3 3 3 3 g 2/4 φ = a cos f 3/27
−
k = 1, 2, 3
−
In the supercritical region when T r
≥ 10, two of these solutions are negative, so
Chen CL
4
the maximal zk is selected as the true compressibility factor. (Note: let z = 0 for C < 0 in the following) Calculate the cpmpressibility factor and molar volume of steam for the reduced temperatures T r = 1 3 and pressures P r = 0.1 10. Prepare a table and a plot of the compressibility factor versus T r and P r as well as a table and a plot of the molar volume versus pressure and T r . The pressure and the volume should be in a logarithmic scale in the second plot.
∼
∼
Chen CL
5
%filename P5_01A_CCL R = 0.08206;%Gas constant (L-atm/g-mol-K) Tc = 647.4; %Critical temperature (K) Pc = 218.3; %Critical pressure (atm) a = 0.42747 * R ^ 2 * Tc ^ (5 / 2) / Pc;%Eq. (4-2), RK equation const b = 0.08664 * R * Tc / Pc; %Eq. (4-3),RK equation consta Pr = 10 %1.2;%0.1; %Reduced pressure (dimensionless) Tr = 3 %1; %Reduced temperature (dimensionless) Asqr = 0.42747 * Pr / (Tr ^ 2.5);%Eq. (4-8) B = 0.08664 * Pr / Tr; %Eq. (4-9) r = Asqr * B; %Eq. (4-6) q = B ^ 2 + B - Asqr; %Eq. (4-7) g = (-27 * r - (9 * q) - 2) / 27;%Eq. (4-12) f = (-3 * q - 1) / 3; %Eq. (4-11) C = (f / 3) ^ 3 + (g / 2) ^ 2; %Eq. (4-10) if (C > 0) E1 = 0 - (g / 2) - sqrt(C); %Eq. (4-15) D = (0 - (g / 2) + sqrt(C)) ^ (1 / 3);%Eq. (4-14) E = sign(E1) * abs(E1) ^ (1 / 3);%Eq. (4-15) z = D + E + 1 / 3;%Eq. (4-13), Compressibility factor (dimensionl else z = 0;
Chen CL
end P = Pr * Pc; T = Tr * Tc; V = z * R * T / P Pr = 10 Tr = 3 V = 0.0837
6
%Pressure (atm) %Temperature (K) %Molar volume (L/g-mol)
Chen CL
%filename P5_01 %filename P5_01B_CCL B_CCL clea cl ear, r, cl clc, c, fo form rmat at co comp mpac act, t, fo form rmat at sh shor ort t g Tc = 64 647. 7.4; 4; %C %Cri riti tica cal l te temp mper erat atur ure e (K (K) ) Pc = 21 218. 8.3; 3; %C %Cri riti tica cal l pr pres essu sure re (a (atm tm) ) Tr_s Tr _set et=[ =[1 1 1. 1.2 2 1. 1.5 5 2 3] 3]; ; Pr_s Pr _set et(1 (1) ) = 0. 0.1; 1; Pr_s Pr _set et(2 (2) ) = 0. 0.2; 2; i = 2; while Pr_se Pr_set(i)< t(i)<=10 =10 i=i+1; Pr_set(i)=Pr_set(i-1)+0.2; end n_Tr n_T r = siz size(T e(Tr_s r_set, et,2); 2); n_Pr n_P r = siz size(P e(Pr_s r_set, et,2); 2); for i=1 i=1:n_ :n_Tr Tr Tr=Tr_set(i); for j=1 j=1:n_ :n_Pr Pr Pr=Pr_set(j); [z(j,i [z( j,i), ), V(j V(j,i) ,i)] ] = RKf RKfun( un(Tr, Tr,Pr, Pr,Tc, Tc,Pc) Pc) ; if z(j z(j,i) ,i)==0 ==0 disp di sp([ ([’ ’ No solu soluti tion on obta obtain ined ed for for Tr = ’ an and Pr Pr = ’ nu num2str(Pr)]);
7
’ nu num m2s 2str tr( (Tr Tr) ) .. ... .
Chen CL
8
end end end disp di sp(’ (’ Co Comp mpre ress ssib ibil ilit ity y Fa Fact ctor or Ve Vers rsus us Tr an and d Pr Pr’) ’); ; disp(’ dis p(’ Tab Tabula ular r Res Result ults’) s’); ; disp(’ dis p(’ ’); disp(’Pr\Tr 1. 1.0 1.2 1.5 2 3 ’); Res=[Pr_se Res=[ Pr_set’ t’ z]; disp(Res); subplot(1,2,1) plot(Pr_set,z(:,1),’-’,Pr_set,z(:,2),’+’,Pr_set,z(:,3),’*’,... Pr_set,z(:,4),’x’,Pr_set,z(:,5),’o’,’LineWidth’,2); set(gca,’FontSize’,14,’Linewidth’,2) legend(’Tr=1’,’Tr=1.2’,’Tr=1.5’,’Tr=2’,’Tr=3’); title( tit le(’\b ’\bf f Com Compre pressi ssibil bility ity Fac Factor tor Ver Versus sus Tr and Pr’ Pr’,’F ,’Font ontSiz Size’, e’,12) 12) xlabel(’\b xlabe l(’\bf f Reduc Reduced ed Press Pressure ure (Pr)’ (Pr)’,’Fon ,’FontSize tSize’,14) ’,14); ; ylabel(’\b ylabe l(’\bf f Compr Compressib essibility ility Facto Factor r (z)’, (z)’,’Font ’FontSize’ Size’,14); ,14); dis di sp( p(’ ’ Pa Paus use e; Plea Pl ease se pres press s an any y ke key y to co cont ntin inue ue ... ... ’) pause P_set=Pr_set.*Pc; disp di sp(’ (’ Mo Mola lar r Vo Volu lume me Ve Vers rsus us Tr an and d P’ P’); ); disp(’ dis p(’ Tab Tabula ular r Res Result ults’) s’); ;
Chen CL
disp(’ ’); disp(’ disp(’Pr\Tr 1. 1.0 1.2 1.5 2 3 ’); Res=[P_set Res=[ P_set’ ’ V]; disp(Res); subplot(1,2,2) loglog(P_set,V(:,1),’-’,P_set,V(:,2),’+’,P_set,V(:,3),’*’,... P_set,V(:,4),’x’,P_set,V(:,5),’o’,’LineWidth’,2); legend(’Tr=1’,’Tr=1.2’,’Tr=1.5’,’Tr=2’,’Tr=3’); set(gca,’FontSize’,14,’Linewidth’,2) title( tit le(’\b ’\bf f Mol Molar ar Vol Volume ume Ver Versus sus Tr and P’, P’,’Fo ’FontS ntSize ize’,1 ’,12) 2) xlabel(’\b xlabe l(’\bf f Press Pressure ure (atm) (atm)’,’Fo ’,’FontSiz ntSize’,14 e’,14); ); ylabel(’\b ylabe l(’\bf f Molar Volum Volume e (L/g(L/g-mol)’ mol)’,’Fon ,’FontSize tSize’,14) ’,14); ;
9
Chen CL
10
function [z, V] = RKfun(Tr,Pr,Tc,Pc) R = 0.08206; %Gas constant (L-atm/g-m a = 0.42747 * R ^ 2 * Tc ^ (5 / 2) / Pc;%Eq. (4-2), RK equation const b = 0.08664 * R * Tc / Pc; %Eq. (4-3),RK equation constant Asqr = 0.42747 * Pr / (Tr ^ 2.5); %Eq. (4-8) B = 0.08664 * Pr / Tr; %Eq. (4-9) r = Asqr * B; %Eq. (4-6) q = B ^ 2 + B - Asqr; %Eq. (4-7) g = (-27 * r - (9 * q) - 2) / 27; %Eq. (4-12) f = (-3 * q - 1) / 3; %Eq. (4-11) C = (f / 3) ^ 3 + (g / 2) ^ 2; %Eq. (4-10) if (C > 0) E1 = 0 - (g / 2) - sqrt(C); %Eq. (4-15) D = (0 - (g / 2) + sqrt(C)) ^ (1 / 3); %Eq. (4-14) E = sign(E1) * abs(E1) ^ (1 / 3); %Eq. (4-15) z = D + E + 1 / 3; %Eq. (4-13), Compressibility factor (dimens else z = 0; end P = Pr * Pc; %Pressure (atm) T = Tr * Tc; %Temperature (K) V = z * R * T / P; %Molar volume (L/g-mol)
Chen CL
11
Calculation of Flow Rate In A Pipeline Concepts Utilized: Application of the general mechanical energy balance for
incompressible fluids, and calculation of flow rate in a pipeline for various pipe diameters and lengths. Numerical Methods: Solution of a single nonlinear algebraic equation and
alternative solution using the successive substitution method. Problem Statement:
The following figure shows a pipeline that delivers water at a constant temperature T = 60oF from point 1 where the pressure is P 1 = 150 psig and the elevation is z1 = 0 ft to point 2 where the pressure is atmospheric and the elevation is z2 = 300 ft.
The density and viscosity of the water can be calculated from the following
Chen CL
12
equations. ρ = 62.122 + 0.0122T 1.54 1057.51 ln µ = 11.0318 + T + 214.624
−
4
2
× 10− T
+ 2.65
7
3
10
× 10− T − 2.24 × 10−
T 4
−
where T is in oF, ρ is in lbm/ft3 , and µ is in lbm/ft s.
·
(a) Calculate the flow rate q (in gal/min) for a pipeline with effective length of L = 1000 ft and made of nominal 8-inch diameter schedule 40 commercial steel pipe. (Solution: v = 11.61 ft/s, gpm =1811 gal/min) (b) Calculate the flow velocities in ft/s and flow rates in gal/min for pipelines at 60oF with effective lengths of L = 500, 1000, . . . 10, 000 ft and made of nominal 4-, 5-, 6- and 8-inch schedule 40 commercial steel pipe. Use the successive substitution method for solving the equations for the various cases and present the results in tabular form. Prepare plots of flow velocity v versus D and L, and flow rate q versus D and L. (c) Repeat part (a) at temperatures T = 40, 60, and 100oF and display the results in a table showing temperature, density, viscosity, and flow rate.
Chen CL
13
Solution:
The general mechanical energy balance on an incompressible liquid applied to this case yields 1 2 g c∆P f F Lv 2 v + g∆z + +2 = 0 (4 20) 2 ρ D where v is the flow velocity in ft/s, g is the acceleration of gravity given by g = 32.174 ft/s2, ∆z = z 2 z1 is the difference in elevation (ft), gc is a conversion factor (in English units gc = 32.174 ft-lbm/lbf s2), ∆P = P 2 P 1 is the difference in pressure lbm/ft2), f F is the Fanning friction factor, L is the length of the pipe (ft) and D is the inside diameter of the pipe (ft). The use of the successive substitution method requires Equation (4-20) to be solved for v as
−
−
−
v=
g∆z +
·
−
g c∆P ρ
1 2
2
f F Lv D
−
2
The equation for calculation of the Fanning friction factor depends on the Reynold’s number, Re = vρD/µ, where µ is the viscosity in lbm/ft-s. For laminar flow (Re < 2100), the Fanning friction factor can be calculated from the equation f F =
16 Re
Chen CL
14
For turbulent flow (Re > 2100) the Shacham equation can be used
1 ε/D f F = log 16 3.7
−
5.02 ε/D 14.5 log + Re 3.7 Re
2
where ε/D is the surface roughness of the pipe (ε = 0.00015 ft for commercial steel pipes). The flow velocity in the pipeline can be converted to flow rate by multiplying it by the cross section area of the pipe, the density of water ( 7.481 gal/ft3), and factor (60 s/min). Thus q has units of (gal/min). The inside diameters (D) of nominal 4-, 5-, 6-, and 8-inch schedule 40 commercial steel pipes are provided in Appendix “Table D-5”. The problem is set up first for solving for one length and one diameter value. The iteration function of the successive substitution method for calculation of the flow velocity is given by vi+l = F (vi), i = 0, l , . . . An error estimate at iteration i is provided by εi = vi
| − vi | +1
The solution is acceptable when the error is small enough, typically εi < 1
5
× 10− .
Chen CL
15
function P5_2C_CCL clear, clc, format short g, format compact D_list=[4.026/12 5.047/12 6.065/12 7.981/12]; % Inside diameter of pi T = 60; %Temperature (deg. F) for i = 1:4 D = D_list(i); j=0; for L=500:500:10000 j = j+1; L_list(j)=L; % Effective length of pipe (ft) [v(j,i),fval]=fzero(@NLEfun,[1 20],[],D,L,T); if abs(fval)>1e-10 disp([’No Conv. for L = ’ num2str(L) ’ and D = ’ num2str(D)]); end q(j,i) = v(j,i) * pi * D ^ 2 / 4* 7.481 * 60; %Flow rate (gpm) end end disp(’ Flow Velocity (ft/s) versus Pipe Length and Diameter’); disp(’ Tabular Results’); disp(’’); disp(’ L\D D=4" D=5" D=6" D=8"’); Res=[L_list’ v];
Chen CL
disp(Res); subplot(1,2,1) plot(L_list,v(:,1),’-’,L_list,v(:,2),’+’,L_list,v(:,3),’*’,... L_list,v(:,4),’x’,’LineWidth’,2); set(gca,’FontSize’,14,’Linewidth’,2) legend(’ D=4"’,’ D=5"’,’ D=6"’,’ D=8"’); title(’\bf Flow Velocity’,’FontSize’,12) xlabel(’\bf Pipe Length (ft)’,’FontSize’,14); ylabel(’\bf Velocity (ft/s)’,’FontSize’,14); disp(’ Pause; Please press any key to continue ... ’) pause disp(’ Flow Rate (gpm) versus Pipe Length and Diameter’); disp(’ Tabular Results’); disp(’’); disp(’ L\D D=4" D=5" D=6" D=8"’); Res=[L_list’ q(:,1) q(:,2) q(:,3) q(:,4)]; disp(Res); subplot(1,2,2) plot(L_list,q(:,1),’-’,L_list,q(:,2),’+’,L_list,q(:,3),’*’,... L_list,q(:,4),’x’,’Linewidth’,2); set(gca,’FontSize’,14,’Linewidth’,2) legend(’ D=4"’,’ D=5"’,’ D=6"’,’ D=8"’); title(’\bf Flow rate’,’FontSize’,14)
16
Chen CL
17
xlabel(’\bf Pipe Length (ft)’,’FontSize’,14); ylabel(’\bf Flow rate (gpm)’,’FontSize’,14);
function fv = NLEfun(v,D,L,T) epsilon = 0.00015;%Surface rougness of the pipe (ft) rho = 62.122+T*(0.0122+T*(-0.000154+T*(0.000000265-... (T*0.000000000224)))); %Fluid density (lb/cu. ft deltaz = 300; %Elevation difference (ft) deltaP = -150; %Pressure difference (psi) vis = exp(-11.0318 + 1057.51 / (T + 214.624)); %Fluid viscosity (lbm/ pi = 3.1416; %The constant pi eoD = epsilon / D; %Pipe roughness to diameter ratio (dimensionless) Re = D * v * rho / vis; %Reynolds number (dimesionless) if (Re < 2100) %Fanning friction factor (dimensionless) fF = 16 / Re; else fF=1/(16*log10(eoD/3.7-(5.02*log10(eoD/3.7+14.5/Re)/Re))^2); end fv=v-sqrt((32.174*deltaz+deltaP*144*32.174/rho)... /(0.5-(2*fF*L/D))); %velocity (ft/s)
Chen CL
18
L\D 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500
Flow Velocity (ft/s) versus Pipe Length and Diameter Tabular Results D=4" D=5" D=6" D=8" 10.773 12.516 14.15 17.035 7.4207 8.6048 9.7032 11.613 5.9721 6.9243 7.8051 9.3295 5.1188 5.9361 6.6912 7.9953 4.5409 5.2674 5.9382 7.0953 4.1168 4.7769 5.3861 6.4362 3.7888 4.3975 4.9592 5.927 3.5255 4.093 4.6166 5.5185 3.3082 3.8416 4.3338 5.1815 3.1249 3.6297 4.0953 4.8973 2.9677 3.4478 3.8907 4.6535 2.8309 3.2896 3.7128 4.4415 2.7106 3.1504 3.5561 4.2548 2.6036 3.0266 3.4169 4.0889 2.5077 2.9156 3.292 3.9402 2.4211 2.8154 3.1793 3.8059 2.3424 2.7244 3.0769 3.6838 2.2706 2.6412 2.9832 3.5723 2.2046 2.5648 2.8972 3.4698
Chen CL
19
10000 2.1437 2.4943 2.8179 3.3752 Pause; Please press any key to continue ... Flow Rate (gpm) versus Pipe Length and Diameter Tabular Results L\D D=4" D=5" D=6" D=8" 500 427.5 780.53 1274.2 2656.4 1000 294.46 536.59 873.81 1811 1500 236.98 431.8 702.87 1454.8 2000 203.12 370.17 602.56 1246.8 2500 180.19 328.48 534.75 1106.4 3000 163.36 297.89 485.03 1003.7 3500 150.35 274.23 446.59 924.25 4000 139.9 255.24 415.74 860.55 4500 131.27 239.56 390.27 807.99 5000 124 226.35 368.8 763.68 5500 117.76 215.01 350.37 725.66 6000 112.33 205.14 334.35 692.6 6500 107.56 196.46 320.24 663.49 7000 103.31 188.74 307.7 637.62 7500 99.508 181.82 296.46 614.43 8000 96.073 175.57 286.31 593.48 8500 92.951 169.89 277.08 574.45
Chen CL
20
9000 9500 10000
90.099 87.48 85.064
164.71 159.94 155.54
268.65 260.91 253.76
557.05 541.07 526.33
Chen CL
21
Adiabatic Operation of A Tabular Reactor for Cracking of Acetone Concepts Utilized: Calculation of the conversion and temperature profile in an
adiabatic tubular reactor. Demonstration of the effect of pressure and heat capacity change on the conversion in the reactor. Numerical Methods:
Solution of simultaneous ordinary differential equations. Problem Statement:
The irreversible, vapor-phase cracking of acetone (A) to ketene (B) and methane (C) that is given by the reaction CH 3COCH 3
→ CH CO + CH 2
4
is carried out adiabatically in a tubular reactor. The reaction is first order with respect to acetone and the specific reaction rate can be expressed by ln(k) = 34.34
− 34222 T
(4-26)
Chen CL
22
where k is in s−1 and T is in K . The acetone feed flow rate to the reactor is 8000 kg/hr, the inlet temperature is T = 1150 K and the reactor operates at the constant pressure of P = 162 kPa (1.6 atm). The volume of the reactor is 4 m3. The material balance equations for the plug-flow reactor are given by dF A = dV dF B = dV dF C = dV
rA (4-27)
−rA (4-28) −rA (4-29)
where F A, F B , and F C are flow rates of acetone, ketene, and methane in g-mol/s, respectively and rA is the reaction rate of A in g-mol/m3 s. The reaction is first order with respect to acetone, thus
·
rA =
−kC A
where C A is the concentration of acetone in g-mol/m 3. For a gas phase reactor, using the appropriate units of the gas constant, the concentration of the acetone in g-mol/m3 is obtained by 1000yAP C A = 8.31T
Chen CL
23
The mole fractions of the various components are given by F A F B F C yA = , yB = , yC = F A + F B + F C F A + F B + F C F A + F B + F C The conversion of acetone can be calculated from
xA =
F A0 F A F A0
−
An enthalpy (energy) balance on a differential volume of the reactor yields dT rA( ∆H ) = dV F AC pA + F B C pB + F C C pC
− −
where ∆H is the heat of the reaction at temperature T (in J/g-mol) and C pA, CpB, and C pC are the molar heat capacities of acetone, ketene and methane (in J/gmol K). Fogler provides the following equations for calculating the heat of
·
Chen CL
24
reaction and the molar heat capacities. ∆H = 80770 + 6.8(T 298) 0.00575(T 2 C pA = 26.60 + 0.183T 45.86 10−6T 2 C pB = 20.04 + .0945T 30.95 10−6T 2 C pC =
2
6
3
3
− − − 298 ) − 1.27 × 10− (T − 298 ) − × − × 13.39 + 0.077T − 18.71 × 10− T 6
2
Galculate the final conversion and the final temperature of P = 1.6, 1.8, . . . , 5.0 atm, for acetone feed rates of F A0 = 10, 20, 30, 35, and 38.3 mol/s where nitrogen is fed to maintain the total feed rate 38.3 mol/s in all cases. Present the results in tabular form and prepare plots of final conversion versus P and F A0 and final temperature versus P and F A0. function P5_03C_CCL clear, clc, format short g, format compact FA0set = [10 20 30 35 38.3]; %Feed rate of acetone in kg-mol/s P_set(1)=1.6; i=1; while P_set(i)<=5; i=i+1; P_set(i)=P_set(i-1)+0.2; % Pressure in atm
Chen CL
25
end n_P = size(P_set’); for i=1:5 FA0=FA0set(i); for j=1:n_P P=P_set(j)*101.325; % Pressure in kPa y0=[FA0; 0; 0; 1035; 0]; [V,y] = ode45(@ODEfun,[0 4],y0,[],FA0,P); Xfin(j,i)=y(end,5); Tfin(j,i)=y(end,4); end end % ------------------------------------------------------------------disp(’ Final Conversion versus FA0 and Pressure’); disp(’ Tabular Results’); disp(’’); disp(’ Pressure FA0=10 FA0=20 FA0=30 FA0=35 FA0=38.3 ’); Res=[P_set’ Xfin(:,1) Xfin(:,2) Xfin(:,3) Xfin(:,4) Xfin(:,5)]; disp(Res); subplot(1,2,1) plot(P_set,Xfin(:,1),’-’,P_set,Xfin(:,2),’+’,P_set,Xfin(:,3),’*’,... P_set,Xfin(:,4),’x’,P_set,Xfin(:,5),’o’,’LineWidth’,2);
Chen CL
26
set(gca,’FontSize’,14,’Linewidth’,2) legend(’FA0=10’,’FA0=20’,’FA0=30’,’FA0=35’,’FA0=38.3’); title(’\bf Final Conversion versus FA0 and Pressure’,’FontSize’,12) xlabel(’\bf Pressure (atm)’,’FontSize’,14); ylabel(’\bf Final Conversion’,’FontSize’,14); disp(’ Pause; Please press any key to continue ... ’) pause disp(’ Final Temperature versus FA0 and Pressure’); disp(’ Tabular Results’); disp(’’); disp(’ Pressure FA0=10 FA0=20 FA0=30 FA0=35 FA0=38.3 ’); Res=[P_set’ Tfin(:,1) Tfin(:,2) Tfin(:,3) Tfin(:,4) Tfin(:,5)]; disp(Res); subplot(1,2,2) plot(P_set,Tfin(:,1),’-’,P_set,Tfin(:,2),’+’,P_set,Tfin(:,3),’*’,... P_set,Tfin(:,4),’x’,P_set,Tfin(:,5),’o’,’Linewidth’,2); set(gca,’FontSize’,14,’Linewidth’,2) legend(’FA0=10’,’FA0=20’,’FA0=30’,’FA0=35’,’FA0=38.3’); title(’\bf Final Temperature versus FA0 and Pressure’,’FontSize’,12) xlabel(’\bf Pressure (atm)’,’FontSize’,14); ylabel(’\bf Temperature (K)’,’FontSize’,14); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chen CL
27
function dYfuncvecdV = ODEfun(V,Yfuncvec,FA0,P); FA = Yfuncvec(1); FB = Yfuncvec(2); FC = Yfuncvec(3); T = Yfuncvec(4); XA = Yfuncvec(5); k = 8.2E14 * exp(-34222 / T); %Reaction rate constant in s-1 FN = 38.3 - FA0; %Feed rate of nitrogene in kg-mol/s yA = FA / (FA + FB + FC + FN); %Mole fraction of acetone CA = yA * P * 1000 / (8.31 * T); %Concentration of acetone in k-mol/m yB = FB / (FA + FB + FC + FN); %Mole fraction of ketene yC = FC / (FA + FB + FC + FN); %Mole fraction of methane rA = -k * CA; %Reaction rate in kg-mole/m3-s deltaH = 80770 + 6.8 * (T - 298) - .00575 * (T ^ 2 - 298 ^ 2)... - 1.27e-6 * (T ^ 3 - 298 ^ 3); CpA = 26.6 + .183 * T - 45.86e-6 * T ^ 2; %Heat capacity of acetone CpB = 20.04 + 0.0945 * T - 30.95e-6 * T ^ 2; %Heat capacity of ketene CpC = 13.39 + 0.077 * T - 18.71e-6 * T ^ 2; %Heat capacity of methane CpN = 6.25 + 8.78e-3 * T - 2.1e-8 * T ^ 2; %Heat capacity of nitrogen dFAdV = rA; %Differential mass balance on acetone dFBdV = -rA; %Differential mass balance on ketene dFCdV = -rA; %Differential mass balance on methane
Chen CL
28
dTdV = (-deltaH)*(-rA)/(FA*CpA+FB*CpB+FC*CpC+FN*CpN); %Differential enthalpy balance dXAdV = -rA / FA0; %Conversion of acetone dYfuncvecdV = [dFAdV; dFBdV; dFCdV; dTdV; dXAdV];
Pressure 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2
Final Conversion versus FA0 and Pressure Tabular Results FA0=10 FA0=20 FA0=30 FA0=35 0.31358 0.27759 0.26394 0.25964 0.32043 0.28346 0.26946 0.26505 0.32651 0.28867 0.27436 0.26986 0.33197 0.29335 0.27877 0.27419 0.33693 0.2976 0.28277 0.27811 0.34147 0.30149 0.28643 0.2817 0.34565 0.30508 0.2898 0.28501 0.34952 0.3084 0.29293 0.28808 0.35313 0.31149 0.29584 0.29094 0.3565 0.31438 0.29856 0.29361 0.35967 0.3171 0.30112 0.29612 0.36266 0.31966 0.30353 0.29849 0.36548 0.32208 0.30581 0.30073 0.36815 0.32438 0.30797 0.30285
Chen CL
4.4 4.6 4.8 5 Column 6 0.25729 0.26265 0.26741 0.27169 0.27558 0.27913 0.28241 0.28545 0.28828 0.29092 0.29341 0.29575 0.29796 0.30006 0.30206 0.30396 0.30578
29
0.37069 0.37311 0.37543 0.37764
0.32656 0.32863 0.33062 0.33251
0.31002 0.31198 0.31385 0.31563
0.30486 0.30678 0.30862 0.31037
Chen CL
30
0.30752 Pause; Please press any key to continue ... Final Temperature versus FA0 and Pressure Tabular Results Pressure FA0=10 FA0=20 FA0=30 FA0=35 Columns 1 through 5 1.6 911.84 908.14 907.47 907.47 1.8 909.02 905.32 904.66 904.67 2 906.5 902.81 902.16 902.18 2.2 904.24 900.55 899.91 899.93 2.4 902.18 898.49 897.87 897.89 2.6 900.3 896.61 895.99 896.02 2.8 898.56 894.87 894.26 894.29 3 896.94 893.26 892.66 892.69 3.2 895.44 891.76 891.16 891.2 3.4 894.03 890.35 889.76 889.8 3.6 892.7 889.03 888.44 888.48 3.8 891.45 887.78 887.19 887.24 4 890.27 886.59 886.02 886.07 4.2 889.15 885.47 884.9 884.95 4.4 888.08 884.4 883.84 883.89 4.6 887.07 883.39 882.82 882.88
Chen CL
31
4.8 5
886.09 885.16
882.42 881.49
881.86 880.93
881.92 880.99
Chen CL
Column 6 907.53 904.74 902.25 900.01 897.97 896.1 894.38 892.78 891.29 889.89 888.58 887.34 886.17 885.05 883.99 882.98 882.02 881.1
32
Chen CL
33
Correlation of The Physical Properties of Ethane Concepts Utilized:
Correlations for heat capacity, vapor pressure, and liquid viscosity for an ideal gas. Numerical Methods:
Polynomial, multiple linear, and nonlinear regression of data with linearization and transformation functions. Problem Statement:
“Tables F-l” through “F-4” of ‘Appendix F” present values for different properties of ethane (ideal gas heat capacity, vapor pressure, and liquid viscosity) as function of temperature. Various regression models will be fitted to the properties of Appendix using MATLAB. (a) Construct a MATLAB function which solves the linear regression problem Xb = y , where X is the matrix of the independent variable values, y is the vector of dependent variable values, and b is the vector of the linear regression model parameters. The input parameters of the function are X , y , and a logical variable which indicates whether there is a free parameter. The returned parameters are: β and the respective confidence intervals, the calculated values of the dependent variable y calc, the linear correlation coefficient R2, and the
Chen CL
34
variance. Test the function by fitting the Wagner equation to vapor pressure data of ethane from “TableF-3” of “Appendix F”. (b) Fit 3rd- and 5th-degree polynomials to the heat capacity of ethane for for the data given in “Tables F-l” and “F-2” of Appendix F by using the multiple linear regression function developed in (a). Compare the quality of the representation of the various data sets with the polynomials of different degrees. (c) Fit the Antoine equation to liquid viscosity of ethane given in “Table F-4” of Appendix F.
Chen CL
35
% filename P5_04A_CCL clear, clc, format short g, format compact prob_title = ([’Vapor Pressure Correlation for Ethane’]); ind_var_name=[’\bf Functions of Reduced Temp.’]; dep_var_name=[’\bf Logarithm of Reduced Pressure’]; fname=input(’Please enter the data file name > ’); % ’VPfile.txt’ xyData=load(fname); X=xyData(:,2:end); y=xyData(:,1); [m,n]=size(X); freeparm=input(’Input 1 if there is a free par., otherwise input 0>’) [Beta,ConfInt, ycal,Var, R2]=MlinReg(X,y,freeparm); disp([’ Results,’ prob_title]); Res=[]; if freeparm==0, nparm = n-1; else nparm = n; end for i=0:nparm if freeparm==0; ii=i+1; else ii=i; end Res=[Res; ii Beta(i+1) ConfInt(i+1)]; end disp(’ Parameter No. Beta Conf_int’); disp(Res); disp([’ Variance ’, num2str(Var)]);
Chen CL
36
disp([’ Correlation Coefficient ’, num2str(R2)]); disp(’ Pause; Please press any key to continue ... ’) pause subplot(1,2,1) plot(y,y-ycal,’*’,’Linewidth’,2) % residual plot set(gca,’FontSize’,14,’Linewidth’,2) title([’\bf Residual, ’ prob_title],’FontSize’,12) xlabel([dep_var_name ’\bf (Measured)’],’FontSize’,14) ylabel(’\bf Residual’,’FontSize’,14) disp(’ Pause; Please press any key to continue ... ’) pause subplot(1,2,2) plot(X,ycal, ’r-’,X,y,’bo’,’Linewidth’,2) title([’\bf Cal/Exp Data ’ prob_title],’FontSize’,12) set(gca,’FontSize’,14,’Linewidth’,2) xlabel([ind_var_name],’FontSize’,14) ylabel([dep_var_name],’FontSize’,14) % VPfile.txt
:
data file provided elsewhere
Please enter the data file name
>
’VPfile.txt’
Chen CL
Input 1 if there is a free parameter, otherwise input 0> Results,Vapor Pressure Correlation for Ethane Parameter No. Beta Conf_int 1 -6.4585 0.09508 2 1.2895 0.21514 3 -1.6712 0.26773 4 -1.2599 0.29417 Variance 9.3486e-005 Correlation Coefficient 1 Pause; Please press any key to continue ... Pause; Please press any key to continue ...
37
0
Chen CL
38
% filename P5_04B_CCL clear, clc, format short g, format compact prob_title = ([’Heat Capacity of Ethane’]); ind_var_name=[’\bf Normalized Temp’]; dep_var_name=[’\bf Heat Capacity J/kmol*K ’]; fname=input(’Please enter the data file name > ’); % ’CPfile.txt’ xyData=load(fname); X=xyData(:,2:end); y=xyData(:,1); [m,n]=size(X); freeparm=input(’Input 1 if there is a free par., otherwise input 0>’) [Beta, ConfInt,ycal, Var, R2]=MlinReg(X,y,freeparm); disp([’ Results, ’ prob_title]); Res=[]; if freeparm==0, nparm = n-1; else nparm = n; end for i=0:nparm if freeparm, ii=i+1; else ii=i; end Res=[Res; ii Beta(i+1) ConfInt(i+1)]; end disp(’ Parameter No. Beta Conf_int’); disp(Res); disp([’ Variance ’, num2str(Var)]);
Chen CL
39
disp([’ Correlation Coefficient ’, num2str(R2)]); plot(y,y-ycal,’*’,’Linewidth’,2) set(gca,’FontSize’,14,’Linewidth’,2) title([’\bf Residual, ’ prob_title],’FontSize’,12) % residual plot xlabel([dep_var_name ’\bf (Measured)’],’FontSize’,14) ylabel(’\bf Residual’,’FontSize’,14) disp(’ Pause; Please press any key to continue ... ’) pause plot(X(:,2),ycal, ’r-’,X(:,2),y,’bo’,’Linewidth’,2) set(gca,’FontSize’,14,’Linewidth’,2) title([’\bf Cal/Exp Data ’ prob_title],’FontSize’,12) xlabel([ind_var_name],’FontSize’,14) ylabel([dep_var_name],’FontSize’,14) Please enter the data file name > ’Cpfile.txt’ Input 1 if there is a free parameter, 0 otherwise > Results, Heat Capacity of Ethane Parameter No. Beta Conf_int 1 34267 1176.6 2 -48870 22701 3 1.0826e+006 1.3785e+005 4 -2.1962e+006 3.4361e+005
1
Chen CL
5 1.8628e+006 3.7147e+005 6 -5.8886e+005 1.4465e+005 Variance 121322.8071 Correlation Coefficient 0.99995 Pause; Please press any key to continue ...
40
Chen CL
41
Complex Chemical Equilibrium by Gibbs Energy Minimization Concepts Utilized:
Formulation of the chemical equilibrium problem as a Gibbs energy minimization problem with atomic balance constraints. The use of Lagrange multipliers to introduce the constraints into the objective function. Conversion of the minimization problem into a system of nonlinear algebraic equations. Numerical Methods:
Solution of a system of nonlinear algebraic equations. Problem Statement:
Ethane reacts with steam to form hydrogen over a cracking catalyst at a temperature of T = 1000 K and pressure of P = 1 atm. The feed contains 4 moles of H 2O per mole of CH 4. Balzisher et al. suggest that only the compounds shown in Table 4-10 are present in the equilibrium mixture (assuming that no carbon is deposited). The Gibbs energies of formation of the various compounds at the temperature of the reaction (1000K) are also given in Table 4-10. The equilibrium composition of the effluent mixture is to be calculated using these data.
Chen CL
42
Table 4-10: Compounds Present in Effluent of Steam Cracking Reactor No.
Comp.
Gibbs Energy kcal/gm-mol
Feed gm-mol
Effluent Ini. Est.
1
CH 4
4.61
0.001
2
C 2H 4
28.249
0.001
3
C 2H 2
40.604
0.001
4
CO2
0.993
5
CO
6
O2
−94.61 −47.942
7
H 2
8
H 2O
9
C 2H 6
1.
0.
0.0001
0.
5.992
−46.03 26.13
4
1.
1
0.001
Formulate the problem as a constrained minimization problem. Introduce the constraints into the objective function using Lagrange multipliers and differentiate this function to obtain a system of nonlinear algebraic equations.
Chen CL
43
Solution:
The objective function to be minimized is the total Gibbs energy given by min ni
G = RT
cni
i=1
Goi RT
+ ln
ni ni
(4
− 49)
where ni is the number of moles of component i, c is the total number of compounds, R is the gas constant, and Go is the Gibbs energy of pure component i at temperature T . The minimization of Equation (4-49) must be carried out subject to atomic balance constraints Oxygen Balance 0 = g 1 = 2n4 + n5 + 2n6 + n7
−4
(4
Hydrogen Balance 0 = g 2 = 4n1 + 4n2 + 2n3 + 2n7 + 2n8 + 6n9 14 Carbon Balance 0 = g 3 = n1 + 2n2 + 2n3 + n4 + n5 + 2n9 2
−
−
− 50) (4 − 51) (4 − 52)
The identification of the various components is given in Table 4-10. These three constraints can be introduced into the objective functions using Lagrange multipliers: λ1, λ2,λ3. The extended objective function is min F =
ni,λj
cni
i=1
3
Goi ni + ln + λj gj RT ni j =1
(4
− 53)
Chen CL
44
The condition for minimum of this function at a particular point is that all the partial derivatives of F with respect to ni and λj vanish at this point. The partial derivative of F with respect to n1 for example, is ∂F Go1 n1 = + ln + 4λ2 + λ3 ∂n1 RT ni
(4
− 54)
The other partial derivatives with respect to ni can be obtained similarly. If it is expected that the amount of a particular compound at equilibrium is very close to zero, it is preferable to rewrite the equation in a form that does not require calculation of the logarithm of a very small number. Rearranging Equation (4-54), for example, yields n1
Go1 + 4λ2 + λ3 = 0 niexp RT
−
The partial derivatives of F with respect to λ1, λ2, and λ3 are g1, g2, and g3, respectively.
Chen CL
45
function P5_05A1_CCL clear, clc, format short g, format compact xguess = [10. 10. 10. 5.992 1. 1. 0.993 0.001 ... 0.001 0.001 0.001 0.0001]; % initial guess vector disp(’Variable values at the initial estimate’); fguess = MNLEfun(xguess); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’int2str(i)’ ’num2str(xguess(i))’ ’ num2str(fguess(i))]); end options = optimset(’Diagnostics’,[’off’],’TolFun’,[1e-9],... ’TolX’,[1e-9]); xsolv = fsolve(@MNLEfun,xguess,options); disp(’Variable values at the solution’); fsolv = MNLEfun(xsolv); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’int2str(i)’ ’num2str(xsolv(i))’ ’ num2str(fsolv(i))]) end
Chen CL
46
function fx = MNLEfun(x) lamda1 = x(1); lamda2 = x(2); lamda3 = x(3); H2 = x(4); H2O = x(5); CO = x(6); CO2 = x(7); CH4 = x(8); C2H6 = x(9); C2H4 = x(10); C2H2 = x(11); O2 = x(12); R = 1.9872; sum = H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2; fx(1,1) = 2 * CO2 + CO + 2 * O2 + H2O - 4; %Oxygen balance fx(2,1) = 4*CH4+4*C2H4+2*C2H2+2*H2+2*H2O+6*C2H6-14; %Hydrogen balance fx(3,1) = CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6 - 2; %Carbo fx(4,1) = log(H2 / sum) + 2 * lamda2; fx(5,1) = -46.03 / R + log(H2O / sum) + lamda1 + 2 * lamda2; fx(6,1) = -47.942 / R + log(CO / sum) + lamda1 + lamda3; fx(7,1) = -94.61 / R + log(CO2 / sum) + 2 * lamda1 + lamda3;
Chen CL
47
fx(8,1) = 4.61 / R + log(CH4 / sum) + 4 * lamda2 + lamda3; fx(9,1) = 26.13 / R + log(C2H6 / sum) + 6 * lamda2 + 2 * lamda3; fx(10,1) = 28.249 / R + log(C2H4 / sum) + 4 * lamda2 + 2 * lamda3; fx(11,1) = C2H2-exp(-(40.604/R+2*lamda2+2*lamda3))*sum; fx(12,1) = O2 - exp(-2 * lamda1) * sum; Variable values at the initial estimate Variable Value Function Value x1 10 -0.0138 x2 10 0 x3 10 0 x4 5.992 19.5944 x5 1 4.6407 x6 1 -6.3214 x7 0.993 -19.8127 x8 0.001 43.2161 x9 0.001 84.0454 x10 0.001 65.1117 x11 0.001 0.001 x12 0.0001 9.9981e-005 Optimization terminated: no further progress can be made. Trust-region radius less than 2*eps.
Chen CL
48
Problem may be ill-conditioned or Jacobian may be inaccurate. Try using exact Jacobian or check Jacobian for errors. Variable values at the solution Variable Value Function Value x1 10+7.44066e-015i -0.013798+4.9091ex2 10-3.7436e-008i -0.015082-0.0002223 x3 10-1.4705e-008i -0.0054569-7.4103ex4 5.992+8.1803e-010i 19.5947+4.04813e x5 1+8.1803e-010i 4.6411+4.0488e-006i x6 1+8.182e-010i -6.3211+4.109e-006i x7 0.993+8.182e-010i -19.8124+4.10898e x8 0.00029165+6.9845e-010i 41.9842+6.3 x9 -0.00037496-3.7053e-005i 83.0696-3. x10 -5.5221e-009-1.0635e-009i 53.0235x11 0.0010003+8.1796e-010i 0.0010003+8 x12 0.00010027+8.1823e-010i 0.00010025
Chen CL
49
function P5_05A2_CCL clear, clc, format short g, format compact xguess = [10. 10. 10. 5.992 1. 1. 0.993 0.001... 0.001 0.001 0.001 0.0001]; % initial guess vector disp(’Variable values at the initial estimate’); fguess = MNLEfun(xguess); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’int2str(i)’ ’num2str(xguess(i))’ ’ num2str(fguess(i))]); end options = optimset(’Diagnostics’,’off’,’TolFun’,1e-9,... ’TolX’,1e-16,’NonlEqnAlgorithm’,’gn’); xsolv = fsolve(@MNLEfun,xguess,options); disp(’Variable values at the solution’); fsolv=MNLEfun(real(xsolv)); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’ int2str(i)’ ’num2str(real(xsolv(i)))’ ’ num2str(fsolv(i)) end
Chen CL
50
function fx = MNLEfun(x) lamda1 = x(1); lamda2 = x(2); lamda3 = x(3); H2 = x(4); H2O = x(5); CO = x(6); CO2 = x(7); CH4 = x(8); C2H6 = x(9); C2H4 = x(10); C2H2 = x(11); O2 = x(12); R = 1.9872; sum = H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2; fx(1,1) = 2 * CO2 + CO + 2 * O2 + H2O - 4; fx(2,1) =4*CH4+4*C2H4+2*C2H2+2*H2+2*H2O+6*C2H6-14; %Hydrogen balance fx(3,1) = CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6 - 2; %Carbo fx(4,1) = log(H2 / sum) + 2 * lamda2; fx(5,1) = -46.03 / R + log(H2O / sum) + lamda1 + 2 * lamda2; fx(6,1) = -47.942 / R + log(CO / sum) + lamda1 + lamda3; fx(7,1) = -94.61 / R + log(CO2 / sum) + 2 * lamda1 + lamda3;
Chen CL
51
fx(8,1) = 4.61 / R + log(CH4 / sum) + 4 * lamda2 + lamda3; fx(9,1) = 26.13 / R + log(C2H6 / sum) + 6 * lamda2 + 2 * lamda3; fx(10,1) = 28.249 / R + log(C2H4 / sum) + 4 * lamda2 + 2 * lamda3; fx(11,1) = C2H2-exp(-(40.604/R+2*lamda2+2*lamda3))*sum; fx(12,1) = O2 - exp(-2 * lamda1) * sum; Variable values at the initial estimate Variable Value Function Value x1 10 -0.0138 x2 10 0 x3 10 0 x4 5.992 19.5944 x5 1 4.6407 x6 1 -6.3214 x7 0.993 -19.8127 x8 0.001 43.2161 x9 0.001 84.0454 x10 0.001 65.1117 x11 0.001 0.001 x12 0.0001 9.9981e-005 Maximum number of function evaluations exceeded. Increase OPTIONS.Max Variable values at the solution
Chen CL
Variable x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
52
Value
Function Value 24.4197 0 0.25306 1.7764e-015 1.5598 0 5.3452 -1.1102e-016 1.5216 2.1094e-015 1.3885 2.2204e-016 0.54492 -3.3307e-015 0.066564 0 1.6707e-007 1.3323e-015 9.5412e-008 1.3323e-015 3.157e-010 1.4387e-020 7.0058e-021 1.5466e-021
Chen CL
53
function P5_05B_CCL clear, clc, format short g, format compact xguess = [10. 10. 10. 5.992 1. 1. 0.993 0.001... 0.001 0.001 0.001 0.0001]; % initial guess vector disp(’Variable values at the initial estimate’); fguess = MNLEfun(xguess); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’ int2str(i)’ ’ num2str(xguess(i))’ ’ num2str(fguess(i))]) end pote=[0 0 0 2 2 2 2 2 2 2 2 2]; tol=1e-9; maxit=100; derfun=0; print=0; [xsolv,y,dy,info]=conles(@MNLEfun,xguess,pote,[],... tol,maxit,derfun,print); disp(’Variable values at the solution’); fsolv = MNLEfun(xsolv); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’ int2str(i)’ ’num2str(xsolv(i))’ ’ num2str(fsolv(i))])
Chen CL
end H2=xsolv(4); H2O=xsolv(5); CO=xsolv(6); CO2=xsolv(7); CH4=xsolv(8); C2H6=xsolv(9); C2H4=xsolv(10); C2H2=xsolv(11); O2=xsolv(12);
54
Chen CL
55
function fx = MNLEfun(x) lamda1 = x(1); lamda2 = x(2); lamda3 = x(3); H2 = x(4); H2O = x(5); CO = x(6); CO2 = x(7); CH4 = x(8); C2H6 = x(9); C2H4 = x(10); C2H2 = x(11); O2 = x(12); R = 1.9872; sum = H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2; fx(1,1) = 2 * CO2 + CO + 2 * O2 + H2O - 4; fx(2,1) = 4*CH4+4*C2H4+2*C2H2+2*H2+2*H2O+6*C2H6-14; %Hydrogen balance fx(3,1) = CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6 - 2; %Carbo fx(4,1) = log(H2 / sum) + 2 * lamda2; fx(5,1) = -46.03 / R + log(H2O / sum) + lamda1 + 2 * lamda2; fx(6,1) = -47.942 / R + log(CO / sum) + lamda1 + lamda3; fx(7,1) = -94.61 / R + log(CO2 / sum) + 2 * lamda1 + lamda3;
Chen CL
56
fx(8,1) = 4.61 / R + log(CH4 / sum) + 4 * lamda2 + lamda3; fx(9,1) = 26.13 / R + log(C2H6 / sum) + 6 * lamda2 + 2 * lamda3; fx(10,1) = 28.249 / R + log(C2H4 / sum) + 4 * lamda2 + 2 * lamda3; fx(11,1) = C2H2-exp(-(40.604/R+2*lamda2+2*lamda3))*sum; fx(12,1) = O2 - exp(-2 * lamda1) * sum; Variable values at the initial estimate Variable Value Function Value x1 10 -0.0138 x2 10 0 x3 10 0 x4 5.992 19.5944 x5 1 4.6407 x6 1 -6.3214 x7 0.993 -19.8127 x8 0.001 43.2161 x9 0.001 84.0454 x10 0.001 65.1117 x11 0.001 0.001 x12 0.0001 9.9981e-005 Variable values at the solution Variable Value Function Value
Chen CL
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
57
24.4197 0.25306 1.5598 5.3452 1.5216 1.3885 0.54492 0.066564 1.6707e-007 9.5412e-008 3.157e-010 5.4592e-021
0 -1.7764e-015 0 -1.1102e-016 -1.6653e-015 -1.5543e-015 5.5511e-015 -4.4409e-016 -1.6964e-013 -2.589e-013 9.3058e-025 -3.9873e-035
Chen CL
58
function P5_05C_CCL clear, clc, format short g, format compact xguess = [10. 10. 10. 5.992 1. 1. 0.993 0.001... 0.001 0.001 0.001 0.0001]; % initial guess vector disp(’Variable values at the initial estimate’); fguess = MNLEfun(xguess); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’ int2str(i)’ ’num2str(xguess(i))’ ’ num2str(fguess(i))]); end pote=[0 0 0 2 2 2 2 2 2 2 2 2]; tol=1e-9; maxit=100; derfun=0; print=0; [xsolv,y,dy,info] = conles(@MNLEfun,xguess,pote,[],... tol,maxit,derfun,print); disp(’Variable values at the solution’); fsolv = MNLEfun(xsolv); disp(’ Variable Value Function Value’) for i=1:size(xguess,2); disp([’x’ int2str(i)’ ’num2str(xsolv(i))’ ’ num2str(fsolv(i))]); end H2 =xsolv(4);
Chen CL
H2O =xsolv(5); CO =xsolv(6); CO2 =xsolv(7); CH4 =xsolv(8); C2H6=xsolv(9); C2H4=xsolv(10); C2H2=xsolv(11); O2 =xsolv(12); R = 1.9872; sum = H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2; G_O2 = O2 * log(abs(O2 / sum)); G_H2 = H2 * log(H2 / sum); G_H2O = H2O * (-46.03 / R + log(H2O / sum)); G_CO = CO * (-47.942 / R + log(CO / sum)); G_CO2 = CO2 * (-94.61 / R + log(CO2 / sum)); G_CH4 = CH4 * (4.61 / R + log(abs(CH4 / sum))); G_C2H6 = C2H6 * (26.13 / R + log(abs(C2H6 / sum))); G_C2H4 = C2H4 * (28.249 / R + log(abs(C2H4 / sum))); G_C2H2 = C2H2 * (40.604 / R + log(abs(C2H2 / sum))); ObjFun=G_H2+G_H2O+G_CO+G_O2+G_CO2... +G_CH4+G_C2H6+G_C2H4+G_C2H2
59
Chen CL
60
function fx = MNLEfun(x) lamda1 = x(1); lamda2 = x(2); lamda3 = x(3); H2 = x(4); H2O = x(5); CO = x(6); CO2 = x(7); CH4 = x(8); C2H6 = x(9); C2H4 = x(10); C2H2 = x(11); O2 = x(12); R = 1.9872; sum = H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2; fx(1,1) = 2 * CO2 + CO + 2 * O2 + H2O - 4; %Oxygen balance fx(2,1) = 4*CH4+4*C2H4+2*C2H2+2*H2+2*H2O+6*C2H6-14; %Hydrogen balance fx(3,1) = CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6 - 2; %Carbo fx(4,1) = log(H2 / sum) + 2 * lamda2; fx(5,1) = -46.03 / R + log(H2O / sum) + lamda1 + 2 * lamda2; fx(6,1) = -47.942 / R + log(CO / sum) + lamda1 + lamda3; fx(7,1) = -94.61 / R + log(CO2 / sum) + 2 * lamda1 + lamda3;
Chen CL
61
fx(8,1) = 4.61 / R + log(CH4 / sum) + 4 * lamda2 + lamda3; fx(9,1) = 26.13 / R + log(C2H6 / sum) + 6 * lamda2 + 2 * lamda3; fx(10,1) = 28.249 / R + log(C2H4 / sum) + 4 * lamda2 + 2 * lamda3; fx(11,1) = C2H2-exp(-(40.604/R+2*lamda2+2*lamda3))*sum; fx(12,1) = O2 - exp(-2 * lamda1) * sum; Variable values at the initial estimate Variable Value Function Value x1 10 -0.0138 x2 10 0 x3 10 0 x4 5.992 19.5944 x5 1 4.6407 x6 1 -6.3214 x7 0.993 -19.8127 x8 0.001 43.2161 x9 0.001 84.0454 x10 0.001 65.1117 x11 0.001 0.001 x12 0.0001 9.9981e-005 Variable values at the solution Variable Value Function Value
Chen CL
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 ObjFun = -104.34
62
24.4197 0.25306 1.5598 5.3452 1.5216 1.3885 0.54492 0.066564 1.6707e-007 9.5412e-008 3.157e-010 5.4592e-021
0 -1.7764e-015 0 -1.1102e-016 -1.6653e-015 -1.5543e-015 5.5511e-015 -4.4409e-016 -1.6964e-013 -2.589e-013 9.3058e-025 -3.9873e-035