3 views (last 30 days)

Show older comments

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

Torsten
on 27 Mar 2020

Edited: Torsten
on 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.

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

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

Start Hunting!