Kinetic parameter estimation for fixed bed reactor
조회 수: 21 (최근 30일)
이전 댓글 표시
I am trying to find the kinetic parameter for the fixed bed reactor for the trireforming of methane having the catalyst weight of w=.02 gram and the initial molar flow rate n0 and final molar flow rate n_Measur is known (mol/sec). There is three reaction for steam reforming, dry reforming oand water gas shift reaction which i am trying to use with total species are CH4,CO2,H2O,O2,N2 (inert),H2 and CO. The code is divided into 4 parts error function, ode function, main function and parameter function.
%% ---------------- Parameter Estimation Main -----------------------------
format long
par = par_kinetic();
%% load data from experimental (to find k and n)
T= 923.15 % Temperature in kelvin
w_Measu = [0 .02]; % weight of catalyst in gram .02
%% nC_Measu final molar flow rate (mol/min) CH4 CO2 H2O O2 N2 H2 CO
nC_Measu = [3.59134E-06 6.12004E-06 1.30813E-08 0 8.32244E-05 2.06871E-05 9.36085E-06
3.55089E-06 2.33395E-05 1.14685E-06 0 6.36501E-05 1.57435E-05 1.57362E-05
2.38098E-05 1.28571E-05 1.95838E-07 0 4.59249E-05 2.74115E-05 1.94572E-05
2.44912E-05 6.85368E-06 9.35903E-06 0 5.61497E-05 1.54736E-05 2.05413E-05
2.43895E-05 1.35797E-05 2.44423E-06 0 4.28068E-05 3.64595E-05 1.91442E-05
4.39838E-06 1.60087E-05 1.7746E-05 0 7.93945E-05 2.17505E-05 1.48277E-05
3.99639E-05 9.70848E-06 5.51779E-06 0 3.64002E-05 3.71545E-05 2.40885E-05
3.14698E-06 1.47095E-05 9.1016E-06 0 6.82368E-05 1.59402E-05 1.20256E-05
3.82879E-05 1.13429E-05 9.03076E-06 0 3.79525E-05 4.11719E-05 2.30244E-05
2.33371E-05 4.29807E-06 1.20675E-05 0 4.97604E-05 3.41758E-05 1.49129E-05
2.12084E-05 2.26973E-05 5.16778E-06 0 3.25403E-05 3.64627E-05 2.14637E-05
4.15565E-06 1.17225E-05 1.40098E-05 0 8.54965E-05 1.46762E-05 1.73545E-05
3.62277E-07 1.64866E-06 0.000153265 0 5.95614E-06 2.11172E-06 1.0817E-06
4.14227E-05 1.30693E-05 7.91973E-08 0 2.5918E-05 4.07423E-05 2.11796E-05
2.38667E-05 2.41786E-06 1.15062E-05 0 5.91536E-05 2.05798E-05 1.30274E-05
2.24993E-05 1.59107E-05 1.68287E-05 0 4.01769E-05 2.02628E-05 2.09384E-05
2.05472E-05 5.95366E-06 3.22445E-06 0 6.23875E-05 3.47837E-05 1.47167E-05
1.94723E-05 2.08046E-05 3.03531E-07 0 4.21599E-05 3.1577E-05 2.33277E-05
1.87001E-05 1.04099E-05 1.87612E-05 0 4.61311E-05 2.91735E-05 1.76068E-05
];
%% Initial Condition for variable MOLAR FLOW RATE(mol/min) CH4 CO2 H2O O2 N2 H2 CO
n0 = [1.11607E-05 5.58036E-06 1.22768E-05 6.69643E-07 8.19196E-05 0 0
1.11607E-05 2.79018E-05 1.22768E-05 6.69643E-07 5.95982E-05 0 0
2.79018E-05 1.67411E-05 2.23214E-05 2.23214E-07 4.44196E-05 0 0
2.79018E-05 1.67411E-05 2.23214E-06 1.11607E-06 6.36161E-05 0 0
2.79018E-05 1.67411E-05 2.23214E-05 1.11607E-06 4.35268E-05 0 0
1.11607E-05 1.67411E-05 1.22768E-05 2.23214E-07 7.12054E-05 0 0
4.46429E-05 1.67411E-05 1.22768E-05 2.23214E-07 3.77232E-05 0 0
1.11607E-05 1.67411E-05 1.22768E-05 1.11607E-06 7.03125E-05 0 0
4.46429E-05 1.67411E-05 1.22768E-05 1.11607E-06 3.68304E-05 0 0
2.79018E-05 5.58036E-06 2.23214E-05 6.69643E-07 5.51339E-05 0 0
2.79018E-05 2.79018E-05 2.23214E-05 6.69643E-07 3.28125E-05 0 0
1.11607E-05 1.67411E-05 2.23214E-06 6.69643E-07 8.08036E-05 0 0
1.11607E-05 1.67411E-05 2.23214E-05 6.69643E-07 6.07143E-05 0 0
4.46429E-05 1.67411E-05 2.23214E-05 6.69643E-07 2.72321E-05 0 0
2.79018E-05 5.58036E-06 1.22768E-05 2.23214E-07 0.000065625 0 0
2.79018E-05 2.79018E-05 1.22768E-05 2.23214E-07 4.33036E-05 0 0
2.79018E-05 5.58036E-06 1.22768E-05 1.11607E-06 6.47321E-05 0 0
2.79018E-05 2.79018E-05 1.22768E-05 1.11607E-06 4.24107E-05 0 0
2.79018E-05 1.67411E-05 1.22768E-05 6.69643E-07 5.40179E-05 0 0
];
%% Initial Guess for parameter (G1)
p0 = [32
12
118];
p1=p0;
%% Solver - Parameter Optimization
for i = 1 : 10
% minimization step
options = optimset('MaxFunEvals', 5000,'Display', 'iter', 'TolX', 1e-10, 'TolFun', 1e-10, 'MaxIter', 5000);
[pmin, resnorm, exitflag, output] = fminsearch(@(p) er_fun(w_Measu,nC_Measu,n0,p),p0,options);
% options = optimoptions('lsqnonlin','Display','iter','Algorithm','levenberg-marquardt');
% [pmin,resnorm,residual,exitflag,output] = lsqnonlin(@(p) er_fun(w_Measu,nC_Measu,n0,p),p0,[],[],options);
% loop with initial condition
p1 = pmin;
%n_Calc;
i
end
%% display value
disp('Estimated parameters p(i):');
disp(pmin)
disp('residual norm:');
disp(resnorm)
%% Estimated parameter
par.p = pmin;
% %% ode solver
%% Calculate molar flow rate using estimated parameters
wspan = [0 0.02];
% Initialize n_Calc to match the dimensions of n0
n_Calc = zeros(size(n0));
% Loop over each row of n0
for i = 1:size(n0, 1)
[~, n_Calc_temp] = ode45(@(w, n) model_fun(w, n, par.p, par), wspan, n0(i, :)); % dNa/dw =ra
n_Calc(i, :) = n_Calc_temp(end, :); % Store the final concentration
end
%% Display molar flow rate
disp('Molar flow rate (mol/min):');
disp(n_Calc);
The main function is use the parameter function.
%% ----------------------- Model Function ---------------------------------
function dy_Calc = model_fun(t_Calc,y_Calc,p,par)
dy_Calc = zeros(size(y_Calc));
PT = 1;
T=923.15;
R = 8.314;
Ksr=1.198*10^17*exp(-26830/T);
Kdr=6.780*10^18*exp(-31230/T);
Kwgs=10^((2078/T)-2.029);
%% parameter
k1 = p(1);
k2 = p(2);
k3 = p(3);
%FN2=0; % Check this value the inert species
% FT=f(1)+f(2)+f(3)+f(4)+f(5);
FT=sum(y_Calc);
pCH4=y_Calc(1)*PT/FT;
pCO2=y_Calc(2)*PT/FT;
pH2O=y_Calc(3)*PT/FT;
pO2=y_Calc(4)*PT/FT;
pN2=y_Calc(5)*PT/FT;
pH2=y_Calc(6)*PT/FT;
pCO=y_Calc(7)*PT/FT;
%% rate expression
r1=k1*pCH4*(1-((pCO*pH2^3)/(pCH4*pH2O*Ksr)));
r2=k2*pCH4*(1-((pCO^2*pH2^2)/(pCH4*pCO2*Kdr)));
r3=k3*((pCO*pH2O/pH2)-(pCO2/Kwgs));
% dy/dt=ra
dy_Calc(1) = - r1-r2; % CH4
dy_Calc(2) = r3-r2; % CO2
dy_Calc(3) = -r3-r1; % H2O
% dy_Calc(4) = 0; % O2
% dy_Calc(5) = 0; % N2
dy_Calc(6) = 3*r1+2*r2+r3; % H2
dy_Calc(7) = -r3+r1+2*r2; % CO
return
Model function defines the rate of reaction and the ode equation
% ---------------------------- error function ----------------------------
function [er,n_Calc]= er_fun(w_Measu,nC_Measu,n0,p)
par = par_kinetic();
%tspan = t_Measu;
%% Solve model with estimated/chosen parameter
wspan = [0 0.02]; %weight of catalyst=tcalc
% [w_Calc, n_Calc] = ode23s(@(t_Calc,nC_Calc) model_fun(t_Calc,nC_Calc,p,par),wspan,n0);
[w_Calc, n_Calc] = ode45(@(w, n) model_fun(w, n, p, par), wspan, n0); % dNa/dw =ra
%% ncalc i have to giive to equation for ra
er = norm(((n_Calc(:,1) - nC_Measu(1,:)).^2+(n_Calc(:,2) - nC_Measu(2,:)).^2 +(n_Calc(:,3) - nC_Measu(3,:)).^2+(n_Calc(:,6) - nC_Measu(6,:)).^2));
end
error function minimizes the objective.
%% ---------------------parameter value------------------------------------
function par = par_kinetic()
%% Model details
par.nofeqns = 3; % number of model equation
par.nVar = 7;
par.unkPar = 3;
return
However, whenever i am runnning the code for all the calculated molar flow rate its coming the be NaN. I am not able to understand what I am doing wrong. Any help will be greatly appreciated. Thanks !!
댓글 수: 6
Torsten
2024년 2월 13일
편집: Torsten
2024년 2월 13일
ode45 computes the time derivatives first with the initial values for the unknowns and will encounter a division by zero when it computes r3. Shall I run the code and show you that r3 will become NaN ? No, I think you can test it yourself by removing the semicolon behind the line
r3=k3*((pCO*pH2O/pH2)-(pCO2/Kwgs));
And remember that ode45 must exactly supply the output for the flow rates at the times when the measurements were done. Thus setting
wspan = [0 0.02];
for the ode integrator is wrong. wspan must be the 20x1 vector when the measurement data were collected.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Biotech and Pharmaceutical에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!