Temperature ramp then keep it at constant temperature
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hello everyone,
The code I had written was at constant temperature.
I would like to do a temperature ramp then hold it at constant temperature.
For example, 200K to 800K at 5K/min then keep it at 800K for 5 hours. Does anyone know how to do it? Thanks a lot!!
The code I had written is in the following:
clear all
clc
format long
%----parameters----%
global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
a=1;
T0=800;%k Temperature
Rg=8.314;%J/k*mol;
E1=10^5;%J/mol;
E2=10^5;%J/mol;
E_2=10^2;%J/mol;
Eg=10^4;%J/mol
k10=0.1;
k20=0.1;
k2_0=0.1;
kgo=0.1;
Deo=10^-10;
k1 = k10*exp(-E1/(Rg*T0));
k2 = k20*exp(-E2/(Rg*T0));
k_2 = k2_0*exp(-E_2/(Rg*T0));
Ke=(2.01*10^(-6))*T0^2-(1.43*10^(-3))*T0+0.2544;
Kg=kgo*exp(Eg/(Rg*T0));
Cb =0;
Ct=20;
Ce=(7.67131719*10^(-47))*T0^(14.5144484);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=linspace(0,100e-4,5);
t=linspace(0,300,100);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = numel(r);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CC0 = zeros(n,1);
CO0 = ones(n,1);
CD0 = ones(n,1);
CM0 = zeros(n,1);
CA0 = zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CO0(1)=15;
CD0(1)=5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y0 = [CC0;CO0;CD0;CM0;CA0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T, Y] = ode15s(@(t,y)fun(t,y,r,n),t,y0)%T is time
%%%%%%%%PLOT%%%%%%%%%%%%%%%%%%%
plot(t,Y)
set(gca,'YLim',[0 1e-2],'xLim',[0 200]);
xlabel('time')
ylabel('concentration')
function DyDt=fun(t,y,r,n)
global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
CC = zeros(n,1);
CO = zeros(n,1);
CD = zeros(n,1);
CM = zeros(n,1);
CA = zeros(n,1);
DCCDt = zeros(n,1);
DCODt = zeros(n,1);
DCDDt = zeros(n,1);
DCMDt = zeros(n,1);
DCADt = zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DyDt = zeros(5*n,1);
rhalf = zeros(n-1,1);
DCCDr = zeros(n-1,1);
D2CCDr2 = zeros(n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CC = y(1:n);
CO = y(n+1:2*n);
CD = y(2*n+1:3*n);
CM = y(3*n+1:4*n);
CA = y(4*n+1:5*n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rhalf(1:n-1)=(r(1:n-1)+r(2:n))/2;
%%%%%%%B.C.Left%%%%%%%%%%%%%%%%%
DCCDr(1)=0;
D2CCDr2(1) = 1.0/(rhalf(1)-r(1))*(CC(2)-CC(1))/(r(2)-r(1));
%%%%%%%B.C.Right%%%%%%%%%%%%%%%%%
DCCDt(n)=Cb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DCODt(n) =-k1*(CO(n)^2 -( 1/Ke)*CD(n)^2);
DCDDt(n) =k1*(CO(n)^2 - (1/Ke)*CD(n)^2)-k2*CD(n)+k_2*CA(n)*CM(n)*CC(n);
DCMDt(n) =k2*CD(n)-k_2*CA(n)*CM(n)*CC(n);
DCADt(n) =k2*CD(n)-k_2*CA(n)*CM(n)*CC(n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:n-1
DCCDr(i) = ((r(i)-r(i-1))/(r(i+1)-r(i))*(CC(i+1)-CC(i))+(r(i+1)-r(i))/(r(i)-r(i-1))*(r(i)-CC(i-1)))/(r(i+1)-r(i-1));
D2CCDr2(i) = (rhalf(i)*(CC(i+1)-CC(i))/(r(i+1)-r(i))-rhalf(i-1)*(CC(i)-CC(i-1))/(r(i)-r(i-1)))/(rhalf(i)-rhalf(i-1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n-1
De=Deo*exp(a*CM(i)/Ct);
DCCDt(i) =2*De*DCCDr(i)+De*r(i)*D2CCDr2(i)+k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
DCODt(i) =-k1*(CO(i)^2 -( 1/Ke)*CD(i)^2);
DCDDt(i) =k1*(CO(i)^2 - (1/Ke)*CD(i)^2)-k2*CD(i)+k_2*CA(i)*CM(i)*CC(i);
DCMDt(i) =k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
DCADt(i) =k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
end
DyDt = [DCCDt;DCODt;DCDDt;DCMDt;DCADt];
end
0 Commenti
Risposta accettata
Torsten
il 27 Mar 2020
Modificato: Torsten
il 27 Mar 2020
Copy the lines where you define the global parameters to the beginning of function "fun" and remove them in the upper part of the code.
Set the end of time integration to 420 minutes.
Remove the two lines
global ...
Change T0=800 to T0=min(200+5*t,800) if t is in minutes or T0=min(200+t/12,800) if t is in seconds.
That's all.
2 Commenti
Torsten
il 28 Mar 2020
Modificato: Torsten
il 28 Mar 2020
The formula for DCCDr(i) is a generalization of the centered difference formula (u'(r)=(u(r+delta)-u(r-delta))/(2*delta)) for non-equidistamt grids. You copied it incorrectly from my code. Please check again.
For the other part:
t=linspace(...) mustn't be defined in fun, but before the call to ode15s.
I never wrote to include the parameters in the list of input variables to fun. They must all be recalculated there because they depend on T0. Copy the lines from a=1 to Ce=... to the beginning of fun and change T0=800 to T0=min(200+5*t,800).
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Performance and Memory in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!