Temperature ramp then keep it at constant temperature

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
 
  

Answer 1

[Torsten]

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.

[S H]

Hello Torsten, I've followed your instruction. The following was the changes I've made.

function DyDt=fun(t,y,r,n,T0,Rg,E1,E2,E_2,Eg,k10,k20,k2_0,kgo,Kg,Deo,k1,k2,k_2,Ke,Cb,Ct,Ce,a)
%global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
t=linspace(0,420,100);
T0=min(200+5*t,800);

However, it gave me this.

Undefined function or variable 't'.

Error in main (line 10)

T0=min(200+5*t,800);

Do you have any idea about it? Where should I revise?

BTW, the code I've written was followed the example you provided couple years ago. (PDE coupled with ODEs by using method of lines) The discretization schemes you provided in this document https://opus4.kobv.de/opus4-zib/files/505/TR-93-14.pdf , page 9-11 describes the discretization of the first spatial derivative. Because it is written in German, it's difficult for me to read. I was wondering if you know the name of the discretization method in English? Is it called klassischen Approximation? I couldn't find any book in English with this approximation. I was wondering if you could provide further information for me to get a better understanding about this method? Thank you very much!

[Torsten]

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).