Need help with Runge Kutta 4th order code

2 views (last 30 days)
Peter Phung
Peter Phung on 16 Oct 2018
I am trying to plot example 5.3 in my heat transfer textbook using 4th order Runge Kutta and on the graph in the text book, it says that tc happens at 124 seconds. Unfortunately for my code, it shows that tc happens at 123 seconds. If I shift the time by one place, it will make tc = 124 seconds, but then the initial temperature (25 degrees celsius) will be shifted to 1 second instead of 0 seconds which is not what the graph in the textbook depicts. tc is defined as the time that it takes for the temperature to reach 150 seconds. I have attached the PDF of the question and my code is provided below.
clear all; close all; clc
tstart=0;
tend=150;
h = 1; %%change to small step size (0.1 or 0.01) then the result of numerical method will be better
t = tstart:h:tend;
w(1) = 25+273;
t(1) = 0;
T(1) = 25+273;
%fprintf('Step 0: t = %12.8f, w = %12.8f\n', t, w);
for i=2:size(t,2)
k1 = h*f(t(i-1),w(i-1));
k2 = h*f(t(i-1)+h/2, w(i-1)+k1/2);
k3 = h*f(t(i-1)+h/2, w(i-1)+k2/2);
k4 = h*f(t(i-1)+h, w(i-1)+k3);
w(i) = w(i-1) + (k1+2*k2+2*k3+k4)/6;
end
figure (1)
neww = w-273 %subtract 723 kelvins from w matrix to convert kelvins back into degrees
plot(t,neww)
title('Plot of temperature vs. time for example 5.3 (first 150 seconds)')
xlabel('Time (seconds)')
ylabel('Temperature (degrees)')
%%%%%%%%%%%%%%%%%%
function deltaT = f(t,T)
convection = 40
tsurr = 175+273
tinf = 175+273
emissivity = 0.8
const = 5.67*10^-8
roe = 2770
ccc = 875
length = (3*10^-3)/2
deltaT = -(convection*(T-tinf)+(emissivity*const*((T^4)-(tsurr^4))))/(roe*ccc*length);
end

Answers (0)

Community Treasure Hunt

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

Start Hunting!