Hello guys!
For an assignment for chemical reactor engineering we have to solve mole balances to calculate the steady state concentrations in a slurry bubble column reactor. To do this, the Matlab function fsolve can be used. Problem is, with the script and function I wrote, Matlab is not giving me a matrix of concentrations but a single row. My goal was to get for every added meter of column (and hence large reactor size/volume) a row of concentrations (so a matrix).
Do you have any idea what's wrong in my script?
Thanks on forehand! Luc
% Reaction rate constants in [m,L^3/g,Ni . s]
k1 = 1.0*10.^8;k2 = 1.0*10.^8;k3 = 1.0*10.^-1;% Activation energies in [J/mol]
Ea1 = 120000;Ea2 = 120000;Ea3 = 150000;% Adsorption constants in [m,L^3/mol]
KA = 0.32;KH2 = 25;KB = 0.07;KC = 0.01;% Particle loading
Ld = 0.10;% Density of catalyst particle in [g/m,cat^3]
rho_cat = 3000000;% Catalyst efficiency
n = 1;% Hold-ups
Es = 0.1; % Given
Eg = 0.11; % Calculated with Krisha 1991 correlation
% Mass-transfer in [1/s]
kGLaGL = 0.1933; % Calculated with Akita & Yoshida correlation
kLSaLS = 0.01;% Saturation concentration of H2 in [mol/m,L^3 . bar]
m = 3;% Pressure in [Bar]
P = 10;% Gas constant in [J/mol . K]
R = 8.314;% Temperature of reaction in [K]
T = 313;% Volume flow in [m,L^3 / s]
Fv = 0.05;% Diameter column/reactor in [m]
Ar = 1;% Height column/reactor in [m]
Hr = linspace(1,25,25);% Ininital concentration of component A in [mol/m,L^3]
CA0 = 1500;% Initial concentrations in [mol/m,L^3]
C0 = [1500 0 0 0 0 0 0 0];% options for accuracy
options = optimoptions('fsolve','TolFun', 1E-6);% Computing the equations
C = fsolve(@(C)AH2BC(C,k1,k2,k3,Ea1,Ea2,Ea3,KA,KH2,KB,KC,Ld,rho_cat,n,Es,Eg,kGLaGL,kLSaLS,m,P,R,T,Fv,Ar,Hr,CA0),C0,options)function [F] = AH2BC(C,k1,k2,k3,Ea1,Ea2,Ea3,KA,KH2,KB,KC,Ld,rho_cat,n,Es,Eg,kGLaGL,kLSaLS,m,P,R,T,Fv,Ar,Hr,CA0)for x=1:25 % Rate equations
r1 = k1*exp(-Ea1/(R*T))*C(3)*KA*C(1)/(1+KA*C(1)+KH2*C(3)+KB*C(5)+KC*C(7)); r2 = k2*exp(-Ea2/(R*T))*C(3)*KB*C(5)/(1+KA*C(1)+KH2*C(3)+KB*C(5)+KC*C(7)); r3 = k3*exp(-Ea3/(R*T))*KA*(C(1).^2)/(1+KA*C(1)+KH2*C(3)+KB*C(5)+KC*C(7)); % Component A liquid phase
% d(CA*Vl)/dt = d(CA*El*Vr)/dt = d(CA*(1-Eg)*Vr)/dt
% Assuming Vr = Vg + Vl (neglecting volume catalyst/solid)
% Acc = in - out - MT(LS)
F(1) = Fv*CA0 - Fv*C(1) - kLSaLS*Es*Ar*Hr(x)*(C(1)-C(2)); % Component A solid phase
% Acc = MT(LS) - reaction
F(2) = kLSaLS*Es*Ar*Hr(x)*(C(1)-C(2)) - (r1+r3)*Ld*rho_cat*Es*Ar*Hr(x)*n; % Component H2 gas phase
% Assuming no MT(GS) and constant volume flow due to constant pressure
% So balanca can be neglected (weak assumption) CH2,g = CH2,sat = P.m
% d(CH2*Eg*Vr)/dt = dFH2/dz - kGLaGL*Eg*Ar*(P*m - CH2(l)) (ODE)
% Component H2 liquid phase
% Acc = -out + MT(GL) - MT(LS)
F(3) = - Fv*C(3) + kGLaGL*Eg*Ar*Hr(x)*((P*m)-C(3)) - kLSaLS*Es*Ar*Hr(x)*(C(3)-C(4)); % Component H2 solid phase
% Acc = MT(LS) - reaction F(4) = kLSaLS*Es*Ar*Hr(x)*(C(3)-C(4)) - (r1+r2)*Ld*rho_cat*Es*Ar*Hr(x)*n; % Component B liquid phase
% Acc = -out - MT(LS)
F(5) = -Fv*C(5) - kLSaLS*Es*Ar*Hr(x)*(C(5)-C(6)); % Component B solid phase
% Acc = MT(LS) - reaction F(6) = kLSaLS*Es*Ar*Hr(x)*(C(5)-C(6)) - (-r1+r2)*Ld*rho_cat*Es*Ar*Hr(x)*n; % Component C liquid phase
% Acc = -out - MT(LS) F(7) = -Fv*C(7) - kLSaLS*Es*Ar*Hr(x)*(C(7)-C(8)); % Component C solid phase
% Acc = MT(LS) - reaction F(8) = kLSaLS*Es*Ar*Hr(x)*(C(7)-C(8)) - (-r2)*Ld*rho_cat*Es*Ar*Hr(x)*n; end end
Best Answer