Question

Write a function to solve the two-dimensional in Matlab, unsteady heat conduction equation with no internal heat generation on a square domain. The length of each side of the square is 1 m. The function will have the following 4 inputs:

npts                     number of grid points in each coordinate direction

nt                         number of time steps to take

dt                         size of time step (s)

alpha                   thermal diffusivity (m²/s)

Use the following initial and boundary conditions:

Initialize the body to T = 0 ^oC. At time t = 0, boundary conditions are imposed such that the temperatures on the sides of the square vary linearly. Here are the corner temperatures:

bottom left corner:         T =   50 ^oC

bottom right corner:      T = 500 ^oC

top right corner:             T = 350 ^oC

top left corner:               T = 150 ^oC

Write the function definition statement. If you type doc function at the MATLAB command prompt it will give you details on the proper syntax. The syntax is general, showing how to use both inputs and outputs. We will have no outputs, only the four inputs noted above.

Make the name of your function the same as the name of your file. For example, if you call your function heat, then your filename should be heat.m.

Create a variable name for each of the four corner boundary conditions and the initial condition, and set their values.

Noting that the temperature on each edge varies linearly, calculate the temperature at every grid point on the outer boundaries. You’ll need the corner temperatures to do this calculation.

Create an array to contain the “old” temperature field at time level p, and call it told. A convenient way to do this is to use the ones function (look at doc ones), which creates an array of all ones of specified dimensions. An easy way to set the field to the initial condition is to multiply this array by the initial temperature.

If there are any calculated variables that you will use multiple times, calculate them. For example, it’s handy to have a variable npm = npts – 1, where npts is the input number of points in each direction.

Create an array to contain the “current” solution at time level p+1, and call it tnew. An easy way to do initialize it in MATLAB is to simply use

tnew = told;

Create a time loop going from 1 to nt, where nt is the input number of time steps.

Inside the time loop, create two other nested loops, one for the “m” indices (x direction) and one for the “n” indices (y direction). Since we don’t have to re-calculate any values on boundaries, the loops will go from 2 to npm.

For each pair of indices within the loops, calculate the solution at time level p+1 (tnew) as a function of temperatures at time level p (told) using the finite-difference form of the 2-D heat equation.

We are going to do a simple animation of the unsteady results as they evolve. Use the contourf function to create the evolving contour plot. (doc contourf) One of the arguments is the number of contours, and 20 seems to work well for this problem. We want to re-draw the contour plot once per time step, so contourf should be called after the end of the nested spatial loops, but before the end of the time loop. It’s not necessary to use a plot handle as we did in the example in class. Simply call contour each time through the loop.

MATLAB tries to be efficient, so when it sees the contourf function, it will buffer the results rather than plotting them right away. We don’t want it to buffer since we want to create a real-time animation. To prevent buffering, include a drawnow statement on a separate line immediately after the call to contourf. (This is similar to the example we did in class.)

Now that the current time step is complete, we have to copy the new solution to the old one so we can go on to the next iteration, told = tnew. Make this the last line before the end of the time loop.

Add the following to the end of the time loop. The first line adds a “colorbar,” which is a scale showing the correspondence between temperatures and colors. The second line adds a label to the colorbar. Finally, the third line gets rid of the plot axes, since we don’t need them here.

hc = colorbar;

hc.Label.String = ‘Temperature, deg. C’;

axis off

When you’re got your program ready to run, use the following input values:

npts = 40

nt = 1000

dt = 0.1

alpha = 0.001

To run the program, simply type

heat(40, 1000, 0.1, 0.001)

Answer 1

%% Here is the Code of the main function. Just copy it all and paste in your Matlab software, then save it as it is and try for different values of inputs.

function Heat(np,nt,dt,alpha)
n= np ; m= np;
L = 1 ; % Length of plate
dx = L / n ;
Omega = (dt* alpha/(dx^2 ));

% Given Temperatures
T0 = 0 ; % Initial Temperature
T1 = 50 ; % Temperature at bottom left corner
T2 = 500 ; % Temperature at bottom right corner
T3 = 350 ; % Temperature at top right corner
T4 = 150 ; % Temperature at top left corner
T_old = T0*ones(nt,m,n);
% Initial conditions
for k=1:nt
for i = 1:m
for j = 1:n
T(k,i,j) = T_old(k,i,j);
end
end
end

% Corner Conditions
T(1,1,1) = 50;
T(1,m,1) = 500;
T(1,m,n) = 350;
T(1,1,n) = 150;
for k=1:nt
for i=2:m %% Bottom Boundary
  
T(k,i,1) = 50 + 450*(i).*dx./L ;
end
for i=2:m %% top Boundary
  
T(k,i,n) = 150 + 200*(i).*dx./L ;
end
for j=2:n %% left Boundary
  
T(k,1,j) = 50 + 100*(j).*dx./L ;
end
for j=2:m %% Right Boundary
  
T(k,m,j) = 500 - 150*(j).*dx./L ;
end
end

T_old = T;
for k = 1:nt-1;
T_old = T ;
for i = 2:m-1
for j = 2:n-1

T(k+1,i,j)= T_old(k,i,j) + Omega*(T_old(k,i+1,j) - 4*T_old(k,i,j)...
+ T_old(k,i-1,j) +T_old(k,i,j+1) + T_old(k,i,j-1));
%



end
end
  
end
x = linspace(0,L,40);
y=x;
[x,y]=meshgrid(x,y);
Tp=squeeze(T(nt,:,:));
Cs = contourf(y,x,Tp,20,'linewidth',2);
xlabel('Length of wall in m','fontsize',14,'fontweight','bold','fontangle','italic');
ylabel('Height of wall in m','fontsize',14,'fontweight','bold','fontangle','italic');
title('Temperature(\circ C) Distribution in Square plate','fontsize'...
,14,'fontweight','bold','fontangle','italic');
set(gca,'fontsize',14,'fontweight','bold','fontangle','italic');
clabel(Cs,'fontsize',14,'fontweight','bold','fontangle','italic');
end

%%%%%% End of the code %%%%%%%%%%%%%%%%%

%% Here is the result for

npts = 40

nt = 1000

dt = 0.1

alpha = 0.001