How to run the same code over and over with one different detail of an event changing each time (temperature the cooling water turns on).

3 views (last 30 days)
Tom Goodland
Tom Goodland on 17 Jan 2022
Commented: Walter Roberson on 19 Jan 2022
I have to find out the maximum time the cooling water can be turned off before being reinstated to avoid the reactor exploding (happens at T = 500 K or P = 45 atm). I want to use the backbone of the code I used for the previous part of my work and have a function keep on picking temperatures to turn on the cooling water above 455 K. The objective for me is to find out how late I can turn the cooling water on and the reactor not explode. Does anyone know what function or code I could use to have UA turn from 0 (NO COOLING WATER) to 2.77E6 (COOLING WATER ON) running the code again for each different temperature the water is turned on? (represents cooling water in the dT derivative)
I think a loop of some sort would help me run this code over and over again with the UA part of the dT/dt derivative changing at different temperatures. However, i'm very unfamiliar with loops and I would appreciate any help.If you need me to provide the code I've been using for the ode, the derivatives, variables, constants etc. let me know. I've posted below the code I used to make an event for previous work on the same problem where the cooling water turns on at 455K.
Any help would be greatly appreciated, thanks.
function main
tspan = [0 10]; %just as an example
initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
options = odeset('Events',@myEventsFcn);
iflag = 1;
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
iflag = 2;
[t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), [te tspan(end)], ye);
t = [t1;te;t2];
y = [y1;ye;y2];
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
cd = y(:,4);
P = y(:,5);
T_2 = y(:,6);
figure(1)
plot(t,[ca,cb,cs,cd])
figure(2)
plot(t,P)
figure(3)
plot(t,T_2)
end
function [value,isterminal,direction] = myEventsFcn(t,y)
value = y(6) - 455; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction

Answers (1)

Walter Roberson
Walter Roberson on 17 Jan 2022
You already have
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
so you already know how to construct anonymous functions to pass in additional parameters.
You can do the same thing for the target temperature
thermostat_setpoints = 455:2.5:500;
ntemps = length(termostat_setpoints);
for tempidx = 1 : ntemps
current_setpoint = thermostat_setpoints(tempidx);
options = odeset('Events', @(t,y)myEventsFcn(t,y,current_setpoint));
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),1), tspan, initial_conditions, options);
[t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),2), [te tspan(end)], ye);
%and so on
end
function [value,isterminal,direction] = myEventsFcn(t,y,setpoint)
value = y(6) - setpoint; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
I would further point out that if you were to add an additional event that detected "the reactor exploding", especially in combination with the second ode15s call, then instead of doing a for loop, you could do a binary search.
  9 Comments
Walter Roberson
Walter Roberson on 19 Jan 2022
stop_temperature = 500;
thermostat_setpoints = 455:2.5:stop_temperature*(1-eps);
ntemps = length(termostat_setpoints);
for tempidx = 1 : ntemps
current_setpoint = thermostat_setpoints(tempidx);
options = odeset('Events', @(t,y)myEventsFcn(t,y,[stop_temperature,current_setpoint);
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),1), tspan, initial_conditions, options);
if any(ie == 1)
%index 1 --> explosion
break;
end
[t2,y2,te2,ye2,ie2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),2), [te tspan(end)], y1(end,:), options);
if any(ie2 == 1)
%index 1 --> explosion
break;
end
end
function [value,isterminal,direction] = myEventsFcn(t,y,setpoints)
value = y(6) - setpoints; % The value that we want to be zero
isterminal = [1 1]; % Halt integration
direction = [0 +1]; %expode --> always halt; turn on cooling only if temperature is increasing
end

Sign in to comment.

Tags

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!