ODE 45 step size changes the turning point

2 ビュー (過去 30 日間)
Xie Yingjie
Xie Yingjie 2021 年 3 月 6 日
回答済み: Stephan 2021 年 3 月 7 日
I face a problem in ODE 45 that when ever i change my step size which is zspan = [0:0.1:1] in this case, the turning point of the graph changes according to it. If i set my step size at 0.1, the turning point occur at 0.1; step size at 0.5, the turning point occur at 0.5.
It makes me feel like the result that i get is very "rigged" according to the step size indicated.
Is there any way to solve this problem?
function F = Val(z,S)
%A-Fe2O3,B-CH4,C-FeO,D-CO2,E-H2O
Xsf=S(1);
Xgf=S(2);
Tsf=S(3);
Tgf=S(4);
%Explicit Information
%Mass balance parameter in fuel reactor
R=8.314; %Gas constant-J/(mol.K)
d=0.5; %fuel reactor diameter-m
A=(pi*0.25^2); %Cross-sectional area-m2
pg=0.2411; %gas density-kg/m3
vch4=0.0628; %Vol. flow rate of methane-m3/s
Mwch4=16; %Molecular weight of methane-g/mol
Fgof=1.247; %Initial molar flow rate of fuel(methane)-mol/s
pbulk=2200; %Bulk density of reactive material in reactor-kg/m3
ps=2.81; %Molar bulk density of reactive material-mol/m3
kof=1.31*(10^8); %Pre-exponential factor
Ea=219900; %Activation energy-J/mol
Cgof=Fgof/vch4; %Initial concentration of fuel(methane)-mol/m3
Tgof=723; %Initial temperature of fuel(methane)-K
n=0.681; %Order of reaction
Fsof=3; %Initial Fe2O3 feed rate-mol/s
kf=kof*exp(-Ea/(R*Tsf));
Cgf=(Cgof*(1-Xgf)*(Tgof))/((1+(2*Xgf))*Tgf);
%Energy balance parameter of fuel reactor
pparticle=4000; %Density of oxygen carrier-kg/m3
dp=75*(10^-6); %mean particle diameter-m
Ap=pi*(dp^2); %Surface area of a particle-m2
Vp=(pi*(dp^3))/6; %Volume of particle-m3
Sm=(pbulk*Ap)/(pparticle*Vp);
Vbar=vch4/A; %Superficial gas velocity-m/s
mug=(0.001596+((3.439*10^-5)*Tgf)+((-8.14*10^-9)*Tgf^2))/1000; %gas viscosity-kg/(m.s)
Nre=(dp*Vbar*pg)/mug;
kg=(-23.35+(0.1698*Tgf)+((1.893*10^-5)*Tgf^2))*(10^-3); %thermal conductivity of gas-W/(m.K)
h=(0.33*(Nre^(1/3))*kg)/dp;
Fsup=0.5; %Molar flow rate of support material-mol/s
CpA=110.9362+(32.04714*(Tsf/1000))+(-9.192333*(Tsf/1000)^2)+(0.901506*(Tsf/1000)^3)+(5.433677/(Tsf/1000)^2); %Heat capacity of Fe2O3-J/(mol.K)
CpB=-0.703029+(108.4773*(Tgf/1000))+(-42.52157*(Tgf/1000)^2)+(5.862788*(Tgf/1000)^3)+(0.678565/(Tgf/1000)^2); %Heat capacity of CH4-J/(mol.K)
CpC=45.7512+(18.78553*(Tsf/1000))+(-5.952201*(Tsf/1000)^2)+(0.852779*(Tsf/1000)^3)+(-0.081265/(Tsf/1000)^2); %Heat capacity of FeO-J/(mol.K)
CpD=24.99735+(55.18696*(Tgf/1000))+(-33.69137*(Tgf/1000)^2)+(7.948387*(Tgf/1000)^3)+(-0.136638/(Tgf/1000)^2); %Heat capacity of CO2-J/(mol.K)
CpE=30.092+(6.832514*(Tgf/1000))+(6.793435*(Tgf/1000)^2)+(-2.53448*(Tgf/1000)^3)+(0.082139/(Tgf/1000)^2); %Heat capacity of H2O-J/(mol.K)
CpF=146.5551+(35.91295*(Tsf/1000))+(-0.183978*(Tsf/1000)^2)+(0.031409*(Tsf/1000)^3)+(-3.659941/(Tsf/1000)^2); %Heat capacity of Support material-J/(mol.K)
HOR1=360642; %Heat of reaction in fuel reactor-J/mol
%Differential Equation
dXgfdz=(ps*A*kf*(Cgf^n))/(4*Fgof);
dXsfdz=(-ps*A*kf*(Cgf^n))/Fsof;
dTgfdz=(A*(Sm*h*(Tsf-Tgf)))/(Fgof*(CpB+Xgf*(2*CpE+CpD+CpB)));
dTsfdz=(A*(Sm*h*(Tsf-Tgf)+((Fgof/A)*(dXgfdz)*HOR1)))/((Fgof*(CpA+(Xsf*(2*CpC-CpA))))+Fsup*CpF);
F=[dXgfdz; dXsfdz; dTgfdz; dTsfdz];
end
To call out the function file:
XsfO = 0.9;
XgfO = 0;
TsfO = 800;
TgfO = 773;
SO = [XsfO XgfO TsfO TgfO];
zspan = [0:0.1:1];
[z,S] = ode45(@val,zspan,SO);
% Plots
figure(1)
hold on
plot(z,(S(:,1)),'-r','displayname','Xsf')
plot(z,(S(:,2)),'-b','displayname','Xgf')
legend;
ylabel('conversion');
xlabel('Length/m');
hold off
figure(2)
hold on
plot(z,(S(:,3)),'-g','displayname','Tsf')
plot(z,(S(:,4)),'-c','displayname','Tgf')
legend;
ylabel('Temperature(K)');
xlabel('Length/m');
hold off

採用された回答

Stephan
Stephan 2021 年 3 月 7 日
This is not a problem of Matlab, but to the misunderstanding how ode solvers in Matlab work. Read here:
To solve this you should use the following notation:
zspan = [0 1];
instead of
zspan = [0:0.1:1];
since the documentation states:
"...However, the solver does not step precisely to each point specified in tspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points in tspan..."
That means, that if you make the steps to big you get inaccurate turning points, if you make the to small you calculate useless points and waste time on calculations. So you should trust the solver to choose the time steps depending on the problem.

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeThermal Analysis についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by