vector and for-loop issues in ode solver

3 views (last 30 days)
Dear all, i am trying to solve a temperature-dependent kinetic expression in ode solver with concentration problem.
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:
%F Interp to stepsize 5
x=0:5:l
F_fit=interp1(x_F,F,x)
function ddx= odefun(~, x) % x contains C
% x(1:n)=C
%%% How to make C alonz z with ceratin value ?
% Interp F
% C=rho./MW.*3.*(1-F_fit)
for j= 2:length(F_fit) %since C= x(1:input.n)
x(1:n,j)= rho./MW.*3.*(1-F_fit(j));
ddx(1)=MW.*(k0.*(rho).*exp(-Ea./(R.*x(n+1))).*x(1,j)...
-A(1).*(rho).*exp(-Ea./(R.*x(n+1))).*(C_0-x(1,j)));
ddx(2:n-1)=MW.*(k0.*(rho).*exp(-Ea./(R.*x(n+2:2*n-1))).*x(2:n-1,j)...
-A(2:n-1).*(rho).*exp(-Ea./(R.*x(n+2:2*n-1))).*(C_0-x(2:n-1,j)));
ddx(n)= MW.*(k0.*(rho).*exp(-Ea./(R.*x(2*n))).*x(n,j)...
-A(n).*(rho).*exp(-Ea./(R.*x(2*n))).*(C_0-x(n,j)));
end
ddx=ddx';
end
  3 Comments
uki71319
uki71319 on 30 Nov 2022
Dear @Star Strider, for Q1, It needs to exist in the calling script workspace before it is passed as an argument.
Do you mean to put f_eq=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))))) before the function or in odefun ??
But then f_eq has a variable T, which would be simulated in odefun, if i put it before odefun, then the error will say unrecozided T.
And even if I've put f_eq in odefun like this,
function dTdz= odefun(~,T,input,p,f_eq)
p=2
f_eq=1-(exp((0.08-7*10^-4*p).*(160+44.*log(p+3)-(T-273.15)))./(1+exp((0.08-7*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)))))
input.Keq=f_eq./(1-f_eq) % will be a vector of 1*input.n
dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*input.C...
-input.k0./input.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./input.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./input.Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*(input.C_0-input.C));
dTdz=dTdz';
end
then i still get this error:
Unrecognized function or variable 'f_eq'.
Error in Question (line 32)
[z, T] = ode45(@(z,T) odefun(z,T,input,p,f_eq), 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);
Star Strider
Star Strider on 30 Nov 2022
It appears that ‘f_eq’ needs to be an anonymous function.
It and all the arguments it needs will have to be passed to ‘odefun’.

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
  7 Comments
Torsten
Torsten on 2 Dec 2022
Edited: Torsten on 2 Dec 2022
doesn't that mean i need to give up using odesolver? Or is it still using the MOL with ode solver to solve it?
It means that you use your measurement data to solve the energy equation using the MOL approach.
--> I thought this is the way for not using kinetic expression; With kinetic expression term, then i can easily use MOL to solve my equations. Am i wrong? sorry i get confused by your way...
So you know the kinetic parameters for the degassing reaction ? And why then do you need F_degree in your model ? If you have the kinetic parameters, you can calculate F_degree by solving simultaneously the equation for C and T.
My new idea is: since from my experimental data, C is already known for certain points along z based on F_degree, thus i just need to interpolate the C data into C_fit along z with z_stepsize. And thus i think i just need to deal with heat balance eq.
But this is exactly what I suggested: deduce the reaction term in mol/(m^3*s) from your measurement data prior to your simulation and only solve the energy equation with these values fixed. But why do you need K_eq / f_eq for this ?
As far as I see, you have two alternatives:
If you have the kinetic parameters for a general degassing kinetic, you don't need F_degree. You can solve C and T simultaneously and get F_degree from the simulation (assuming that the other model parameters are known and don't need to be determined).
If you don't have the kinetic parameters, then - prior to your simulation - you can estimate the molar fluxes for z (and assuming uniform distribution over r, also for r) from F_degree. Having these terms, you can include them in the energy equation and solve it alone (without using the equation for C).
And with respect to the mentionned aim of your simulation:
I find it quite strange that you want to fit the two parameters you mentionned (thermal conductivity and heat transfer coefficient) by some kind of simulation. There are well-suited experimental designs (especially to measure thermal conductivity).
Torsten
Torsten on 2 Dec 2022
But prior to my simulation --> That means i need to write a for-loop?
I think pencil and paper is more suited regarding the number of elements in F_degree.

Sign in to comment.

More Answers (0)

Categories

Find more on Chemistry in Help Center and File Exchange

Tags

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!