Prev: OPC "trend" function
Next: variables in titles ??
From: Torsten Hennig on 8 Aug 2010 22:44
%%%%%%%%%%%%%% MAIN FUNCTION using Newton Raphson method %%%%%%%%%%%%%%%%
function gasf_Newton_6sJ
global T
global E
global Stoichiometry
global P

% input the flows of elements entering the system
nC = input('enter -> nC, mol = ');
nH = input('enter -> nH, mol = ');
nO = input('enter -> nO, mol = ');
E = [nC; nH; nO];
% set the temperature (K)
T = input('enter -> T, K = ');
% set the pressure (Po = 1 bar)
P = input('enter -> P, bar = ');
%set the stiociometry matrix
Stoichiometry = [1 0 1 0 1 0;...
4 2 0 2 0 0;...
0 1 1 0 2 2];

% solve the system of equations
[N] = Newton(@f,@jac,[1 1 1 1 1 1]',1e-6);

% outputs of system, moles and total moles
CH4 = sqrt(N(1));
H2O = sqrt(N(2));
CO = sqrt(N(3));
H2 = sqrt(N(4));
CO2 = sqrt(N(5));
O2 = sqrt(N(6));
nt = CH4+H2O+CO+H2+CO2+O2;

>> I was wrong in my last reply.
>> The last seven lines should read
>> CH4 = N(1)^2
>> H2O = N(2)^2
>> CO = N(3)^2
>> H2 = N(4)^2
>> CO2 = N(5)^2
>> O2 = N(6)^2
>> nt = CH4+H2O+CO+H2+CO2+O2

% mole fractions
yCH4 = CH4/nt
yH2O = H2O/nt
yCO = CO/nt
yH2 = H2/nt
yCO2 = CO2/nt
yO2 = O2/nt
ntt = yCH4 + yH2O + yCO + yH2 + yCO2 + yO2

%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobian Matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Jac] = jac(N)

sum = N(1)^2+N(2)^2+N(3)^2+N(4)^2+N(5)^2+N(6)^2;

Jac(1:3,:) = [2*N(1), 0, 2*N(3), 0, 2*N(5), 0 ;...
8*N(1), 4*N(2), 0, 4*N(4), 0, 0 ;...
0, 2*N(2), 2*N(3), 0, 4*N(5), 4*N(6)];

% CH4 + H2O -> CO + 3*H2
% lnKp1 = 2*ln(P/nt) + ln(nCO) + 3*ln(nH2) - ln(nCH4) - ln(nH2O)
Jac(4,1) = (-4*N(1))/sum-2/N(1);
Jac(4,2) = (-4*N(2))/sum-2/N(2);
Jac(4,3) = (-4*N(3))/sum+2/N(3);
Jac(4,4) = (-4*N(4))/sum+6/N(4);
Jac(4,5) = (-4*N(5))/sum;
Jac(4,6) = (-4*N(6))/sum;
% CO + H2O -> CO2 + H2
% lnKp2 = ln(nCO2) + ln(nH2) - ln(nCO) - ln(nH2O)
Jac(5,1) = 0;
Jac(5,2) = -2/N(2);
Jac(5,3) = -2/N(3);
Jac(5,4) = 2/N(4);
Jac(5,5) = 2/N(5);
Jac(5,6) = 0;
% H2 + 0.5*O2 -> H2O
% lnKp4 = -0.5*ln(P/nt) + ln(nH2O) - ln(nH2) - 0.5*ln(nO2)
Jac(6,1) = N(1)/sum;
Jac(6,2) = N(2)/sum+2/N(2);
Jac(6,3) = N(3)/sum;
Jac(6,4) = N(4)/sum-2/N(4);
Jac(6,5) = N(5)/sum;
Jac(6,6) = N(6)/sum-1/N(6);

%%%%%%%%%%%%%%%%%%%%% Equilibrium Constants, lnKp %%%%%%%%%%%%%%%%%%%%%%%
function [LnKp1 LnKp2 LnKp3] = LnKp(T)
% database of temperatures, K
Tdb = [300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 ...
1800 1900 2000];

% database of lnKp values from GasEq
lnKp1 = [-56.864 -35.9835 -23.2035 -14.5385 -8.2618 -3.5002 0.2359 3.2434...
5.716 7.7828 9.536 11.04 12.3441 13.4842 14.4895 15.3812 16.1773 16.8919];
lnKp2 = [11.4575 7.3537 4.9302 3.3489 2.2479 1.4448 0.8381 0.3667 -0.0081...
-0.3117 -0.5614 -0.7696 -0.9449 -1.0942 -1.2225 -1.3333 -1.4296 -1.5141];
lnKp3 = [91.6186 67.3317 52.7002 42.9072 35.8863 30.6019 26.479 23.1708...
20.4575 18.191 16.2691 14.6193 13.187 11.9323 10.8234 9.8368 8.9531...
8.1571];

% interpolate equilibrium constant
LnKp1 = interp1(Tdb,lnKp1,T,'spline','extrap');
LnKp2 = interp1(Tdb,lnKp2,T,'spline','extrap');
LnKp3 = interp1(Tdb,lnKp3,T,'spline','extrap');

%%%%%%%%%%%%%%%%%%%%% Set of Equations to Solve %%%%%%%%%%%%%%%%%%%%%%%%%
function [eq] = f(N)
global Stoichiometry
global E
global T
global P

% the element balance equations
eq(1:3) = Stoichiometry*N.^2 - E.^2;

>> If your first version (before transforming
>> the variables) was correct, this should read
>> eq(1:3) = Stoichiometry*N.^2 - E;

% the equilibrium constants
[LnKp1 LnKp2 LnKp3] = LnKp(T);

% equilibrium equations
eq(4) = -LnKp1+2*log(P/sum(N.^2))+log(N(3)^2)+3*log(N(4)^2)-log(N(1)^2)-log(N(2)^2);
eq(5) = -LnKp2+log(N(5)^2)+log(N(4)^2)-log(N(3)^2)-log(N(2)^2);
eq(6) = -LnKp3-0.5*log(P/sum(N.^2))+log(N(2)^2)-log(N(4)^2)-0.5*log(N(6)^2);
eq = eq';

%%%%%%%%%%%%%%%%%%%%%%% Newton Raphson Method %%%%%%%%%%%%%%%%%%%%%%%%%%%
% modified from example in Numercial Methods in Chemical Engineering
function solution = Newton(Function,Jacobian,x_i,tol)
x = x_i;
error = 2*tol;
iter = 0;
while error > tol
F = feval(Function,x);
error1 = max(abs(F));
error2 = max(abs(F));
J = feval(Jacobian,x);
dx = J\(-F);
m = 1;
while error2 >= error1||~isreal(F)
xnew = x+0.5^m*dx;
F = feval(Function,xnew);
error2 = max(abs(F));
m = m+1;
end
x = xnew;
F = feval(Function,x);
error = max(abs(F));
iter = iter+1;
disp(['error = ' num2str(error)])
disp(['iter = ' num2str(iter)])
end
solution = x;

Best wishes
Torsten.
From: Gavin on 9 Aug 2010 05:28
> % outputs of system, moles and total moles
> CH4 = sqrt(N(1));
> H2O = sqrt(N(2));
> CO = sqrt(N(3));
> H2 = sqrt(N(4));
> CO2 = sqrt(N(5));
> O2 = sqrt(N(6));
> nt = CH4+H2O+CO+H2+CO2+O2;
>
> >> I was wrong in my last reply.
> >> The last seven lines should read
> >> CH4 = N(1)^2
> >> H2O = N(2)^2
> >> CO = N(3)^2
> >> H2 = N(4)^2
> >> CO2 = N(5)^2
> >> O2 = N(6)^2
> >> nt = CH4+H2O+CO+H2+CO2+O2
>
> % the element balance equations
> eq(1:3) = Stoichiometry*N.^2 - E.^2;
>
> >> If your first version (before transforming
> >> the variables) was correct, this should read
> >> eq(1:3) = Stoichiometry*N.^2 - E;

Torsten,
I made the changes which seems to have taken care of the imaginary number problem. The model also runs from T = 700K-1500K with no warnings (original model only ran from T = 1200K-1500K with no warnings). However, the model fails to converge at T = 600K, and at T = 300K-500K matlab starts to give warnings again. The warning given at T = 300K is as follows:

> In gasf_Newton_6sJ>Newton at 126
In gasf_Newton_6sJ at 23
error = 1.4951e-06
iter = 162
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.137846e-35.

It would be great if I could avoid these warnings. My goal is to fix this so I can avoid problems when I use the model for a larger system. I've posted the current m-file if you want to run it and see if you get the same warning. My inputs are:
enter -> nC, mol = 2
enter -> nH, mol = 2
enter -> nO, mol = 4
enter -> T, K = 300
enter -> P, bar = 1

latest m-file is here... https://files.me.com/wigging/s2d25v

Gavin
From: Torsten Hennig on 9 Aug 2010 01:55
> > % outputs of system, moles and total moles
> > CH4 = sqrt(N(1));
> > H2O = sqrt(N(2));
> > CO = sqrt(N(3));
> > H2 = sqrt(N(4));
> > CO2 = sqrt(N(5));
> > O2 = sqrt(N(6));
> > nt = CH4+H2O+CO+H2+CO2+O2;
> >
> > >> I was wrong in my last reply.
> > >> The last seven lines should read
> > >> CH4 = N(1)^2
> > >> H2O = N(2)^2
> > >> CO = N(3)^2
> > >> H2 = N(4)^2
> > >> CO2 = N(5)^2
> > >> O2 = N(6)^2
> > >> nt = CH4+H2O+CO+H2+CO2+O2
> >
> > % the element balance equations
> > eq(1:3) = Stoichiometry*N.^2 - E.^2;
> >
> > >> If your first version (before transforming
> > >> the variables) was correct, this should read
> > >> eq(1:3) = Stoichiometry*N.^2 - E;
>
> Torsten,
> I made the changes which seems to have taken care of
> the imaginary number problem. The model also runs
> from T = 700K-1500K with no warnings (original model
> only ran from T = 1200K-1500K with no warnings).
> However, the model fails to converge at T = 600K,
> , and at T = 300K-500K matlab starts to give warnings
> again. The warning given at T = 300K is as follows:
>
> > In gasf_Newton_6sJ>Newton at 126
> In gasf_Newton_6sJ at 23
> error = 1.4951e-06
> iter = 162
> Warning: Matrix is close to singular or badly
> dly scaled.
> Results may be inaccurate. RCOND =
> ccurate. RCOND = 1.137846e-35.
>
> It would be great if I could avoid these warnings.
> My goal is to fix this so I can avoid problems when
> n I use the model for a larger system. I've posted
> the current m-file if you want to run it and see if
> you get the same warning. My inputs are:
> enter -> nC, mol = 2
> enter -> nH, mol = 2
> enter -> nO, mol = 4
> enter -> T, K = 300
> enter -> P, bar = 1
>
> latest m-file is here...
> https://files.me.com/wigging/s2d25v
>
> Gavin

Just out of curiosity:
Why don't you use MATLAB's fsolve for your
nonlinear system of equations ?

Best wishes
Torsten.
From: Gavin on 9 Aug 2010 08:47
I have tried using fsolve but ended up with incorrect results. Feel free to try with my system of equations and let me know if you get it to work correctly. I found that the newton raphson method gave more consistent answers and was easier to fix if problems arose.

Gavin
From: Gavin on 10 Aug 2010 09:57
> > Gavin
>
> Just out of curiosity:
> Why don't you use MATLAB's fsolve for your
> nonlinear system of equations ?
>
> Best wishes
> Torsten.

Torsten,

I created a model using fsolve but still have problems at lower temperature ranges. At higher temperatures it will run but often gives incorrect answers. The warning that it gives at lower temperatures is the following...
>>Optimizer appears to be converging to a point which is not a root.
Norm of relative change in X is less than max(options.TolX^2,eps) but
sum-of-squares of function values is greater than or equal to sqrt(options.TolFun)
Try again with a new starting guess.
I receive much better results when implementing the newton raphson method that we have been discussing.

Does anyone know how I can make the newton raphson model stable at lower temperature ranges?

Gavin
First  |  Prev  | 
Pages: 1 2 3 4 5
Prev: OPC "trend" function
Next: variables in titles ??