vector and for-loop issues in ode solver

2 次查看（过去 30 天）
uki71319
uki71319 2022-11-30
编辑： uki71319 2023-3-27
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 个评论
显示 1更早的评论
uki71319
uki71319 2022-11-30
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 2022-11-30
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’.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-11-30
%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 个评论
显示 5更早的评论
Torsten
Torsten 2022-12-2
编辑：Torsten 2022-12-2
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 2022-12-2
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by