Consider a long insulated rod is initially heated electrically along its length

Question

Consider a long insulated rod is initially heated electrically along its length with both ends in ice bath. This creates the following information:

-boundary condition T=0 when x=0,L for all time =t

-initial condition T=Sin(x) when t=0 for 0<x<L

Now, the electric heater is turned off. Plot the transient temp distribution along the rod as the rod cools off. Use copper as the material with =117x10-6 m2/s (Page 269). L=1m.

Write your computer program (any language you are comfortable, MatLab, Java, Excel, etc.) as a general case to be able to use the program for the following input:

Thermal Diffusivity, ALPHA, Rod Length, L, Initial Tinial=Sin(x), Number of Elements in x-direction ,N.

For saving computational time and memory, consider half the rod, since this is a symmetric problem.

Direction:

Calculate mesh size Dx=L/N and select a dt such that r<0.5 for a stable solution. Eg. For input L=1m and N=10 notes, Dx=0.1 . Also select a dt for a stable solution, i.e. r<0.5 – e.g. dt=0.002

Calculate the initial temp at each node by putting the appropriate value of Dx at each node. i.e. Ti,p= Sin(x) for p=1

March forward in time to get Ti,(p+1) using the node temp data T i,p from above at (p) time. Calculate all the nodes temp: Repeat this process, first by marching in i (x-direction using Dx=0.1) and then march to the next p (time direction using dt=0.002

Plot your data at t=0, .01, .02, .04, .06, .08, .1, .15, .2 and .3

Then consider an unstable case ,i.e. when r>0.5 – Try dt=0.01 and r=1.0

Now plot your data at 0.1, 0.16 and 0.22

Explanation / Answer

i wrote two different types of codes for the temperature of rod varaitaons

electric heater rod temperature matlab code
boundary condition T=0 when x=0,L for all time =t
initial condition T=Sin(x) when t=0 for 0<x<L
*********************************
first code
//
dz = .25; %each depth step is 1/4 meter
Nz = 400; % Choose the number of depth steps (should go to at least 100 m)
Nt = 5000; % Choose the number of time steps
dt = (365*24*60*60)/Nt; %Length of each time step in seconds (~ 6.3*10^3 seconds, or ~105 minutes)
K = 2*10^-6; % "canonical" K is 10e-6 m^2/s

T = 15*ones(Nz+1,Nt+1); %Create temperature matrix with Nz+1 rows, and Nt+1 columns
% Initial guess is that T is 15 everywhere.
time = [0:12/Nt:12];
T(1,:) = 15-10*sin(2*pi*time/12); %Set surface temperature

maxiter = 500
for iter = 1:maxiter
Tlast = T; %Save the last guess
T(:,1) = Tlast(:,end); %Initialize the temp at t=0 to the last temp
for i=2:Nt+1,
depth_2D = (T(1:end-2,i-1)-2*T(2:end-1,i-1)+T(3:end,i-1))/dz^2;
time_1D = K*depth_2D;
T(2:end-1,i) = time_1D*dt + T(2:end-1,i-1);
T(end,i) = T(end-1,i); % Enforce bottom BC
end
err(iter) = max(abs(T(:)-Tlast(:))); %Find difference between last two solutions
if err(iter)<1E-4
break; % Stop if solutions very similar, we have convergence
end
end
if iter==maxiter;
warning('Convergence not reached')
end

figure(1)
plot(log(err)), title('Convergence plot')
figure(2)
imagesc([0 12],[0 100],T); title('Temperature plot (imagesc)')
colorbar
figure(3)
depth = [0:dz:Nz*dz];
contourf(time,-depth,T); title('Temperature plot (contourf)')
colorbar

figure(4)
plot(time,T(1,:),time,T(21,:),time,T(41,:),time,T(61,:),time,T(81,:))
xlabel('Time (months)'); ylabel('Temperature (C)');
legend('0m','5m','10m','15m', '20m')

**********************
second code
//
function ypheatode = ypheatode(t,u)
K = .001; rho = 1.; c = 1.;
L = 1.; dx = L/4; const = (K/(rho*c))/(dx*dx);
f = [1 1 1]/(rho*c);
ypheatode(1) = -const*(2*u(1) - u(2)) + f(1);
ypheatode(2) = -const*(-u(1) + 2*u(2) - u(3)) + f(2);
ypheatode(3) = -const*(2*u(3) - u(2)) + f(3);
ypheatode = [ypheatode(1) ypheatode(2) ypheatode(3)]';

************
clear; clf
K = .001; rho = 1.; c = 1.;
L = 1.; dx = L/4; const = (K/(rho*c))/(dx*dx);
f = [1 1 1]/(rho*c);
uo = [1 1 1];
to = 0;
tf = 400;
[t u] = ode45('ypheatode',[to tf],uo);
maxk = size(t,1);
x = 0:.25:1;
figure(1)
for k = 1:5:maxk
plot(x,[0 u(k,:) 0])
title('temperatures at various times')
xlabel('position in space')
axis([0 1 0 130])
hold on
pause
end
t(1:5:maxk)
A = const*[2 -1 0;-1 2 -1;0 -1 2];
d = [1 1 1]'/(rho*c);
sssol = Ad
plot(x, [0 sssol' 0])
figure(2)
plot(t,u(:,2))
title('temperature in the center')
xlabel('time')
************************************

Related Questions

Navigate

• Browse (All)
• Browse C
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.