vector and for-loop issues in ode solver
显示 更早的评论
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 个评论
Star Strider
2022-11-30
I do not see any assignment to ‘f_eq’ in the posted code.
It needs to exist in the calling script workspace before it is passed as an argument.
uki71319
2022-11-30
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
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 个评论
But what if in the future, i want to use f_eq as anoymous function in ode? How can i modify my code to work?
Why ? I don't see any advantage to do so.
I was thinking that you told me to solve mass balance equation as well, but then the C concentration needs to solved or set at certain point on z direction
Yes, you get a set of discretized equations for C just like the ones you work with for T.
Torsten
2022-12-1
I cannot decipher which PDEs you are solving for T and C.
I described the way how I would determine the amount of developing gas in mol/(m^3*s) for each cell in z- and r- direction. This value can be used in the energy equation for each cell to compute the source term.
I cannot understand your C computations in the new code.
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
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 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Translated by 
