# vector and for-loop issues in ode solver

7 views (last 30 days)
uki71319 on 30 Nov 2022
Commented: uki71319 on 21 Dec 2022
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]
%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);
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);

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
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?

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2022a

### Community Treasure Hunt

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

Start Hunting!

Translated by