Prev: OPC "trend" function
Next: variables in titles ??
From: Gavin on 6 Aug 2010 00:06
> > Gavin
>
> Did you adjust the Jacobian matrix according to the substitution ?
>
> Best wishes
> Torsten.

Torsten,
I have now adjusted the jacobian accordingly. However, I'm having issues with making the rest of the code consistent with changing to the N(i)^2 terms. I've posted links to my files below. Let me know how you would fix this. I am running the model with inputs of nC=2, nH=2, nO=4, T=500..1500, P=1. It runs fine at temperatures greater than 1200K but not so well below 1200K. The answers at 1300K should be yCH4=4.36e-8, yH2O=0.24882, yCO=0.24882, yH2=0.08451, yCO2=0.41785.
Gavin

I've posted my original m-file here... https://files.me.com/wigging/fh1src
and my modified (i.e. N(i)^2 ) file here... https://files.me.com/wigging/47n7xa
From: Torsten Hennig on 6 Aug 2010 01:12
> > > Gavin
> >
> > Did you adjust the Jacobian matrix according to the
> substitution ?
> >
> > Best wishes
> > Torsten.
>
> Torsten,
> I have now adjusted the jacobian accordingly.
> However, I'm having issues with making the rest of
> f the code consistent with changing to the N(i)^2
> terms. I've posted links to my files below. Let me
> know how you would fix this. I am running the model
> with inputs of nC=2, nH=2, nO=4, T=500..1500, P=1.
> It runs fine at temperatures greater than 1200K but
> t not so well below 1200K. The answers at 1300K
> should be yCH4=4.36e-8, yH2O=0.24882, yCO=0.24882,
> yH2=0.08451, yCO2=0.41785.
> Gavin
>
> I've posted my original m-file here...
> https://files.me.com/wigging/fh1src
> and my modified (i.e. N(i)^2 ) file here...
> https://files.me.com/wigging/47n7xa

In every expression in the function F, you will have
to replace N(i) by N(i)^2.
This gives new functions F1,...,F6 which you will
have to partially differentiate with respect to
the N(i) to put the result in function Jac.

Best wishes
Torsten.
From: Gavin on 6 Aug 2010 19:42
> In every expression in the function F, you will have
> to replace N(i) by N(i)^2.
> This gives new functions F1,...,F6 which you will
> have to partially differentiate with respect to
> the N(i) to put the result in function Jac.
>
> Best wishes
> Torsten.

Torsten,
I've made the N(i)^2 changes to the F function but now I'm getting imaginary numbers. Something is definitely not right. I've updated the jacobian accordingly too. I've posted the latest file here:
https://files.me.com/wigging/6f6j7j

Gavin
From: Torsten Hennig on 7 Aug 2010 02:15
I did not check your equations in detail, just
found two errors that I marked.

Best wishes
Torsten.


%%%%%%%%%%%%%% 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 = sqrt(sum(N));

>> nt must be the sum of the sqrt, not the
>> sqrt of the sum.

% 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)

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(N)-2/N(1);
Jac(4,2) = (-4*N(2))/sum(N)-2/N(2);
Jac(4,3) = (-4*N(3))/sum(N)+2/N(3);
Jac(4,4) = (-4*N(4))/sum(N)+6/N(4);
Jac(4,5) = (-4*N(5))/sum(N);
Jac(4,6) = (-4*N(6))/sum(N);
% 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(N);
Jac(6,2) = N(2)/sum(N)+2/N(2);
Jac(6,3) = N(3)/sum(N);
Jac(6,4) = N(4)/sum(N)-2/N(4);
Jac(6,5) = N(5)/sum(N);
Jac(6,6) = N(6)/sum(N)-1/N(6);

>> You still use sum(N) in the Jac-routine although
>> you changed N to N.^2.

%%%%%%%%%%%%%%%%%%%%% 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;
% 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;
From: Gavin on 8 Aug 2010 21:44
> % 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 = sqrt(sum(N));
> >> nt must be the sum of the sqrt, not the
> >> sqrt of the sum.
>
> Jac(6,4) = N(4)/sum(N)-2/N(4);
> Jac(6,5) = N(5)/sum(N);
> Jac(6,6) = N(6)/sum(N)-1/N(6);
> >> You still use sum(N) in the Jac-routine although
> >> you changed N to N.^2.

Torsten,
I've made the changes as you stated above and still get the wrong solutions. Imaginary numbers are still showing up somehow. Uploaded the latest m-file and the jacobian m-file that I'm using.
Gavin

m-file for model... https://files.me.com/wigging/10xkqn
jacobian m-file... https://files.me.com/wigging/cpg7p8
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5
Prev: OPC "trend" function
Next: variables in titles ??