• Author / Uploaded
• Augusto De La Cruz Camayo

Citation preview

Cutlip and Shacham: Problem Solving in Chemical and Biochemical 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 Concepts Utilized: Analytical solution of the cubic Redlich-Kwong equation for 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 RT a √ − (V − b) V (V +!b) T 5/2 R 2 Tc a = 0.42747 Pc   RTc b = 0.08664 Pc

P

=

(1)

Chen CL

2

where

P V T R Te Pe The compressibility

= pressure in atm = molar volume in liters/g-mol = temperature in K = gas constant (R = 0.08206 (atm-liter/g-mol-K)) = critical temperature in K = critical pressure in atm factor is given by PV z= 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 − z 2 − qz − r = 0 where

r = A2B, A2 = 0.42747

q 

Pr 5/2 Tr

 B

(5) = B 2 + B − A2   = 0.0866 PTrr

in which Pr is the reduced pressure (P/Pc) and Tr is the reduced temperature (T /Tc). 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 C f

 3   f g 2 = + 3 2 −3q − 1 −27r − 9q − 2 g = = 3 27

If C > 0, there is one real solution for z given by z = D + E + 1/3 √ 1/3 D = (−g/2 + C)

E = (−g/2 −

√

C)1/3

If C < 0, there are three real solutions r





−f φ 2π(k − 1) 1 cos + + , 3s 3 3 3 g 2/4 φ = a cos −f 3/27

zk

= 2

k = 1, 2, 3

In the supercritical region when Tr ≥ 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 Tr = 1 ∼ 3 and pressures Pr = 0.1 ∼ 10. Prepare a table and a plot of the compressibility factor versus Tr and Pr as well as a table and a plot of the molar volume versus pressure and Tr . 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 consta b = 0.08664 * R * Tc / Pc; %Eq. (4-3),RK equation constan 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 (dimensionle 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_01B_CCL clear, clc, format compact, format short g Tc = 647.4; %Critical temperature (K) Pc = 218.3; %Critical pressure (atm) Tr_set=[1 1.2 1.5 2 3]; Pr_set(1) = 0.1; Pr_set(2) = 0.2; i = 2; while Pr_set(i) 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 (dimensi 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 P1 = 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 × 10−4T 2 + 2.65 × 10−7T 3 − 2.24 × 10−10T 4 1057.51 ln µ = −11.0318 + T + 214.624 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 gc∆P fF Lv 2 1 2 +2 = 0 (4 − 20) − v + g∆z + 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 = z2 − z1 is the difference in elevation (ft), gc is a conversion factor (in English units gc = 32.174 ft-lbm/lbf·s2), ∆P = P2 − P1 is the difference in pressure lbm/ft2), fF 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 s v=

g∆z +

gc∆P ρ

.

Lv 2

1 fF −2 2 D



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 fF =

16 Re

Chen CL

14

For turbulent flow (Re > 2100) the Shacham equation can be used    2 1 ε/D 5.02 ε/D 14.5 fF = log − log + 16 3.7 Re 3.7 Re 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 × 10−5.

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 pip 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/f 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 CH3COCH3 → CH2CO + CH4 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 dFA dV dFB dV dFC dV

=

rA (4-27)

= −rA (4-28) = −rA (4-29)

where FA, FB , and FC 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 = −kCA where CA is the concentration of acetone in g-mol/m3. 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 CA = 8.31T

Chen CL

23

The mole fractions of the various components are given by FB FC FA , yB = , yC = yA = FA + FB + FC FA + FB + FC F A + FB + FC The conversion of acetone can be calculated from

xA =

FA0 − FA FA0

An enthalpy (energy) balance on a differential volume of the reactor yields dT −rA(−∆H) = dV FACpA + FB CpB + FC CpC where ∆H is the heat of the reaction at temperature T (in J/g-mol) and CpA, CpB, and CpC 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 CpA CpB CpC

= = = =

80770 + 6.8(T − 298) − 0.00575(T 2 − 2982) − 1.27 × 10−6(T 3 − 2983) 26.60 + 0.183T − 45.86 × 10−6T 2 20.04 + .0945T − 30.95 × 10−6T 2 13.39 + 0.077T − 18.71 × 10−6T 2

Galculate the final conversion and the final temperature of P = 1.6, 1.8, . . . , 5.0 atm, for acetone feed rates of FA0 = 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 FA0 and final temperature versus P and FA0. 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) ’); % ’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

40

5 1.8628e+006 6 -5.8886e+005

3.7147e+005 1.4465e+005 Variance 121322.8071 Correlation Coefficient 0.99995 Pause; Please press any key to continue ...

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 H2O per mole of CH4. 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

CH4

4.61

0.001

2

C2H4

28.249

0.001

3

C2H2

40.604

0.001

4

CO2

−94.61

0.993

5

CO

−47.942

1.

6

O2

0.

0.0001

7

H2

0.

5.992

8

H2 O

9

C2H6

−46.03

4

1.

26.13

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  o  X Gi G ni min = cni + ln P (4 − 49) ni RT RT ni i=1 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 = g1 = 2n4 + n5 + 2n6 + n7 − 4 (4 − 50) Hydrogen Balance 0 = g2 = 4n1 + 4n2 + 2n3 + 2n7 + 2n8 + 6n9 − 14 (4 − 51) Carbon Balance 0 = g3 = n1 + 2n2 + 2n3 + n4 + n5 + 2n9 − 2 (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

X i=1

 cni

Goi



3

X ni + λj gj + ln P ni RT 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 P + 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 −

X

 niexp

Go1 RT

 + 4λ2 + λ3 = 0

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; %Carbon 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.9091e-0 x2 10-3.7436e-008i -0.015082-0.00022231 x3 10-1.4705e-008i -0.0054569-7.4103e-0 x4 5.992+8.1803e-010i 19.5947+4.04813ex5 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.10898ex8 0.00029165+6.9845e-010i 41.9842+6.35 x9 -0.00037496-3.7053e-005i 83.0696-3.0 x10 -5.5221e-009-1.0635e-009i 53.0235-2 x11 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; %Carbon 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.MaxF 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; %Carbon 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; %Carbon 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

Chen CL

63

Thank You for Your Attention Questions Are Welcome