vector and for-loop issues in ode solver

7 views (last 30 days)
Dear all, i am trying to solve a temperature-dependent kinetic expression in ode solver with concentration problem.
With presumption of C:
Q1: The code with the given temeparure works without putting in ode, but it doesn't work when it is for temperature simulation in ode solver.(Solved)
With the correct way to use C:
Q2: I've decided to solve mass and heat balance eqs. together. Since the concentration C of degas will be depentdent on F_degree along z direction, but ode solver will develop the z automatically --> I try to write a for loop for C with regarding F_degree in ode solver, but the code doesn't work at all....
Can anyone tell me where is the problem and how can i solve it? ? I'm not sure if making a loop for C in z direction(j) in ode is a correct way or not...
Any help or advice is greatly appreciated!!
Codes:
%Parameters
input.n=10; % # in r
input.Ea=63000; % [J/mol]
input.R=8.314; % [J/mol/K]
input.k0=505; % [m^3/kg/s]
input.H=112500; % [J/mol]
input.rho=360; % [kg/m^3]
input.MW=0.3; % [kg/mol]
input.T0=linspace(300,350,input.n)
l=0.8 % [m] z len
p=2 % [bar]
input.r=linspace(0,0.2,input.n) % radius=[m]
%Concentration
%%% Correct calculation %%%
%reaction:
% A-->A_X + 3*X
z_F=linspace(0,0.8,6)
F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85]
C_0=input.rho./input.MW.*3.*(1-F_degree(1)) %[kg/m^3] % C_0: X gas still remained on substnace A at the start
input.C_0=C_0
%%% plot
for i=1:length(F_degree)
C(i)= input.rho./input.MW.*3.*(1-F_degree(i)); % C: X gas still remained on substance A in the process
% remained
end
C
plot(z_F,C,'r.',MarkerSize=10)
xlabel("z [m]"); ylabel("C [mol/m^3]")
title('C vs z')
%F_dgeree Interp to z stepsize 0.005
z=0:0.005:l
F_degree_fit=interp1(z_F,F_degree,z) % with real experimental data
plot(z,F_degree_fit,'b*')
title("z vs F\_degree\_fit)")
% X gas equilibrium Rate
% r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C)
% r = input.k0 *(input.rho)*exp(-input.Ea/(input.R*T))*input.C - input.k0./input.Keq *(input.rho)*exp(-input.Ea/(input.R*T))*(input.C_0-input.C)
% P.S. (C_0-C) means the released X gas
x0= zeros(2.*input.n,1);
x0(1:input.n)=input.C_0 % C_0
x0(input.n+1:2*input.n)=input.T0 % Inlet temperature T0
[z, x] = ode45(@(z,x) odefun(z,x,input), 0:0.005:l, x0)
[R,Z]=meshgrid(input.r,z)
surf(R,Z,x(:,input.n+1:2*input.n))
function ddz= odefun(~, x, input) % x contains C and T
% x(1:input.n)=C
% x(input.n+1:2*input.n)=T
p=2; %[bar]
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))));
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))));
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))));
Keq=f_eq./(1-f_eq) ; % will be a vector of 1*input.n
Keq=Keq.';
%%% How to make C alonz z with ceratin value ????
% F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85];
% Interp F_degree with stepsize: 0.005 --> F_degree_fit
% C=input.rho./input.MW.*3.*(1-F_degree_fit)
for j= 2:length(F_degree_fit) %since C= x(1:input.n)
x(1:input.n,j)= input.rho./input.MW.*3.*(1-F_degree_fit(j));
dCdz(1)=input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1,j)...
-input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1,j)));
dCdz(2:input.n-1)=input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1,j)...
-input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1,j)));
dCdz(input.n)= input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n,j)...
-input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n,j)));
dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1,j)...
-input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1,j)));
dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1,j)...
-input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1,j)));
dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n,j)...
-input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n,j)));
end
ddz=[dCdz,dTdz]';
end
But i got this error:
Unrecognized function or variable 'F_degree_fit'.
Error in Question>odefun (line 64)
for j= 2:length(F_degree_fit) %since C= x(1:input.n)
Error in Question (line 43)
[z, x] = ode45(@(z,x) odefun(z,x,input), 0:0.005:l, x0)
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
  4 Comments
uki71319
uki71319 on 30 Nov 2022
Dear @Star Strider, I then put the anonymous function for f_eq,
[z, T] = ode45(@(z,T) odefun(z,T,input), 0:0.005:l, T0)
function dTdz= odefun(~,T,input)
p=2
f_eq=@(T) 1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T-273.15)))))
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))))
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))))
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))))
Keq=f_eq./(1-f_eq) ; % will be a vector of 1*input.n
Keq=Keq.'
dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*input.C...
-input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*(input.C_0-input.C));
dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*input.C...
-input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*(input.C_0-input.C));
dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*input.C...
-input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*(input.C_0-input.C));
dTdz=dTdz';
end
but it turns out to be :
Conversion to function_handle from double is not possible.
Error in Question>odefun (line 40)
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))))
Error in Question (line 34)
[z, T] = ode45(@(z,T) odefun(z,T,input), 0:0.005:l, T0)
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Sign in to comment.

Accepted Answer

Torsten
Torsten on 30 Nov 2022
%Parameters
input.n=10; % # in r
input.Ea=63000; % [J/mol]
input.R=8.314; % [J/mol/K]
input.k0=505; % [m^3/kg/s]
input.H=112500; % [J/mol]
input.rho=360; % [kg/m^3]
input.MW=0.3; % [kg/mol]
l=0.8 % [m] z total length
p=2 % [bar]
%Concentration
%%% Presumption_1st_try_in_ode, which is wrong
F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85]
% A--> B+3*gas
input.C_0=input.rho./input.MW.*3.*(1-F_degree(1)) % [mol/m^3] the [C] of gas molecule still in A at the start
input.C=input.rho./input.MW.*3.*(1-F_degree(end)) % [mol/m^3] the [C] of gas molecule still in A
%%% Correct calculation for concentration %%%
z=linspace(0,0.8,6)
C_0=input.rho./input.MW.*3.*(1-F_degree(1)) % at the beginning
C=input.rho./input.MW.*3.*(1-F_degree)
%%% plot
for i=1:length(F_degree)
C(i)= input.rho./input.MW.*3.*(1-F_degree(i));
end
C
plot(z,C,'r.',MarkerSize=10)
xlabel("z [m]"); ylabel("C [mol/m^3]")
title('C vs z')
%Rate
%r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C)
%r = input.k0 *(input.rho)*exp(-input.Ea/(input.R*T))*input.C - input.k0./input.Keq *(input.rho)*exp(-input.Ea/(input.R*T))*(input.C_0-input.C)
T0=linspace(300,350,input.n);
[z, T] = ode45(@(z,T) odefun(z,T,input), 0:0.005:l, T0)
function dTdz= odefun(~,T,input)
p=2;
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))));
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))));
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))));
Keq=f_eq./(1-f_eq);
Keq = Keq.';
dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*input.C...
-input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*(input.C_0-input.C));
dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*input.C...
-input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*(input.C_0-input.C));
dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*input.C...
-input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*(input.C_0-input.C));
dTdz=dTdz';
end
  13 Comments
uki71319
uki71319 on 21 Dec 2022
Dear @Torsten, so following your second suggestion,
You know the molar flux entering (ndot_in) and leaving (ndot_out) the reactor and the residence time T in the reactor.
Thus (ndot_in - ndot_out)*T [moles] of your substance degas.
You know the free volume V of the reactor.
Thus (ndot_in - ndot_out)*T/V [moles/m^3] degas in total
--> since i could directly use the experimental C [mol/m^3] from degree of F along z, so these sentesnces (ndot_in - ndot_out)*T/V [moles/m^3] turns out to be : CA_in - CA_out, [mol/m^3] right?
But then,
Now using your F, you can distribute this total amount (ndot_in - ndot_out)*T/V to the different z-zones of the reactor.
For z-Zone i : (ndot_in - ndot_out)*T/V * deltaFi.
where I assume that the deltaFi sum to 1 (maybe in your case, you have to set ndot_out = 0).
Now you can distribute this total amount for z-zone i to the r-zones in z-Zone i:
For r-Zone j in z-Zone i : (ndot_in - ndot_out)*T/V * deltaFi * pi*(r_(j+1)^2-r_j^2) / (pi*R^2)
-> How come these above make it to the final reaction unit: [mol/(m^3*s)] ? Since F_degree is no unit.
Now multiply this expression with -delta_r H and you'll arrive at an approximate value for the heat released in cell (i,j)
--> and later if I keep going to multiply (CA_in - CA_out) with (-delta_r H),
then the final unit become: [mol/(m^3)]*[J/mol] = [J/m^3] , instead of what we want: [J/(m^3*s)]...
Thank you for making these detail instructions before, but when i follow it throgh, i found it still confusing. Could you please tell me more about it?

Sign in to comment.

More Answers (0)

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by