Step change for ODE

Hi,

I have a function that solves the dynamic behavior of a reactor over time given certain input conditions using the ODE solver. However, I would like to be able to see what would happen if a step change in the input condition or some other variable is applied at a certain t. I believe Simulink can be used to determine the response to a step function, but I'm not sure how to pass my function to Simulink. For example, in my code below, the input condition is Tc. I would like to see the change if the variable q was changed from 150 to 100 at t = 20 if Tc is the same or if Tc was changed at t = 20.

Thanks for any help!

function ivp
% Inputs
tval=0:0.01:20;
Tc = [290; 300; 305; 310];
% Differential Equation
Ca0=0.5; T0=350; Tc0 = 300;
Int = [Ca0; T0];
[t,Y1] = ode15s(@ode,tval,Int,[],Tc(1));
% Plots
figure()
plot(t,Y1(:,2),'r')
xlabel('Time (min)'); ylabel('Reactor temperature (K)');legend('290K');
figure()
plot(t,Y1(:,1),'r')
xlabel('Time (min)');ylabel('Reactant A concentration (mol/L)');legend('290K');
end

function [Dy] = ode(t,y,Tc)
q=150; Cai=1; Ti=350; V=100; p=1000; C=0.239; H=5*10^4;
ER=8750; k0=7.2*10^10; UA=5*10^4; 
% y(1) is Ca   Dy(1) is deriv
% y(2) is T    Dy(2) is deriv
Dy  = zeros(2,1);
Dy(1) = q/V*(Cai-y(1))-k0*exp(-ER/y(2))*y(1);
Dy(2) = q/V*(Ti-y(2))+H/(p*C)*k0*exp(-ER/y(2))*y(1)+UA/(V*p*C)*(Tc-y(2));
end

Answer 1

[Star Strider]

Try this:

function ivp
% Inputs
tval=0:0.01:20;
Tc = [290; 300; 305; 310];
% Differential Equation
Ca0=0.5; T0=350; Tc0 = 300;
Int = [Ca0; T0];
[t1,Y11] = ode15s(@ode,tval,Int,[],Tc(1));
tval = tval+20.01;                                              % New ‘t’
Int = Y11(end,:);                                               % Last Values Are New Initial Conditions
[t2,Y12] = ode15s(@ode,tval,Int,[],Tc(2));                      % New Value For ‘Tc’
t = [t1; t2];                                                   % Concatenate Time Vectors
Y1 = [Y11; Y12];                                                % Concatenate Integrated Outputs
% Plots
figure()
plot(t,Y1(:,2),'r')
xlabel('Time (min)'); ylabel('Reactor temperature (K)');legend('290K');
figure()
plot(t,Y1(:,1),'r')
xlabel('Time (min)');ylabel('Reactant A concentration (mol/L)');legend('290K');
end

function [Dy] = ode(t,y,Tc)
% q=150; 
q = 150.*(t<=20) + 100.*(t>20);                                 % Define ‘q’ With Step Change At ‘t=20’
Cai=1; Ti=350; V=100; p=1000; C=0.239; H=5*10^4;
q = 150.*(t<=20) + 100.*(t>20);
ER=8750; k0=7.2*10^10; UA=5*10^4; 
% y(1) is Ca   Dy(1) is deriv
% y(2) is T    Dy(2) is deriv
Dy  = zeros(2,1);
Dy(1) = q/V*(Cai-y(1))-k0*exp(-ER/y(2))*y(1);
Dy(2) = q/V*(Ti-y(2))+H/(p*C)*k0*exp(-ER/y(2))*y(1)+UA/(V*p*C)*(Tc-y(2));
end

This changes ‘Tc’ in the second ode15s call, and changes ‘q’ inside the ‘ode’ function. Experiment with other options.

EDIT — (9 Aug 2019 at 21:59)

To plot the temperatures and their effects on the plot, with corresponding legend descriptions, the code changes:

function ivp
% Inputs
tval=0:0.01:20;
Tc = [290; 300; 305; 310];
% Differential Equation
Ca0=0.5; T0=350; Tc0 = 300;
Int = [Ca0; T0];
[t1,Y11] = ode15s(@ode,tval,Int,[],Tc(1));
tval = tval+20.01;                                              % New ‘t’
Int = Y11(end,:);                                               % Last Values Are New Initial Conditions
[t2,Y12] = ode15s(@ode,tval,Int,[],Tc(2));                      % New Value For ‘Tc’
% t = [t1; t2];                                                   % Concatenate Time Vectors
% Y1 = [Y11; Y12];                                                % Concatenate Integrated Outputs
% Plots
figure()
plot(t1,Y11(:,2),'r')
hold on
plot(t2,Y12(:,2),'g')
hold off
xlabel('Time (min)')
ylabel('Reactor temperature (K)')
legend('290K','300K');
figure()
plot(t1,Y11(:,1),'r')
hold on
plot(t2,Y12(:,1),'g')
hold off
xlabel('Time (min)')
ylabel('Reactant A concentration (mol/L)')
legend('290K','300K');
end

If you have more temperatures and more values for ‘q’, a loop will be most efficient, and a different definition for ‘q’ will be necessary. This code works for two temperatures and two values for ‘q’.

Except for the plot section changes, my code is unchanged fro m my original Answer.

[Dennis T]

Thanks so much!

[Star Strider]

As always, my pleasure!

Answer 2

[Fahad Jetpurwala]

Good evening. Can someone please help asap? Ive got a model for a pem fuel cell in which im having troubles understanding the correct way to put in step input for 4 of my variables and then finding the step response of one of the input variables. This is the code for matlab. Please explain what im doing wrong? Thanks in advance. time=[1 30 60];%time period F_air=[ 100 250 500];%%air flow rate 500ml/min step(time,F_air) title('Step Response of Air Flow rate'); figure F_Fuel=[ 500 800 1000];%%fuel flow rate%1000ml/min step(time,F_Fuel) title('Step Response of Fuel Flow rate'); figure P_Fuel=[ 0.02 0.2 0.4];%%fuel Pressure MPA 0.02-0.4 step(time,P_Fuel) title('Step Response of Fuel Pressure'); figure P_air=[ 0.2 1 2];%%Air Pressure MPA 0.2-2 step(time,P_air) title('Step Response of Air Pressure');

%Constant Values of PEM Fuel Cell

FuelFlowRate= 1000;%fuel flow rate%1000ml/min AirFlowRate = 500;%%air flow rate 500ml/min FuelPressure= 0.45;%fuel Pressure MPA 0.02-0.4 AirPressure = 2;%Air Pressure MPA 0.2-2

%% output plots figure plot(current) figure plot(voltage) figure plot(stack) figure plot(H2FlowRate)

[Star Strider]

Please post this as a new Question, and format it correctly using the Code radio button in the top toolstrip.

It does not belong here.

I will delete it from here tomorrow.