LAB 1: MATLAB - Introduction to Programming

LAB 1: MATLAB - Introduction to Programming
Objective:
The objective of this laboratory is to review how to use MATLAB as a programming tool and to
review a classic analytical solution to a steady-state two-dimensional heat transfer conduction
problem. We will also review how to write a program that is well documented and easy to read.
Background:
The steady-state two-dimensional temperature distribution, T(x, y), in a solid with constant
thermal conductivity and rectangular cross section as shown in Figure 1 can be found
analytically by solving the heat diffusion equation (or Laplace’s Equation)
∂ 2T ∂ 2T
+
=0
∂x 2 ∂y 2
(1)
using the Separation of Variables method. An introduction to the method can be found in
Bergman, et al., Introduction to Heat Transfer, 6th ed., Wiley, 2011, pp 231-235. Another
€ solutions to conduction problems is Carslaw and Jaeger,
excellent reference for analytical
Conduction of Heat in Solids, Oxford University Press, 1959. We are interested in analytical
solutions in this class because they are useful for validating a numerical method (or making sure
our approximate numerical calculations are accurate).
y
y = Ly
T(x, y = Ly) = T2
T(x = 0, y) = T1
T(x = Lx, y) = T1
x
T(x, y = 0) = T1
x = Lx
Figure 1. Schematic of two-dimensional domain for conduction heat transfer.
As an example, for the following boundary conditions shown in Figure 1
T ( x = 0, y ) = T1, T ( x = Lx , y ) = T1, T ( x, y = 0) = T1, T ( x, y = Ly ) = T2
(2)
where T1 and T2 are constant temperatures and dimensionless temperature θ ( x, y ) is defined as
€
€
€
€
T ( x, y ) − T1
θ ( x, y ) =
T2 − T1
(3)
€
the analytical solution can be determined to be
€
2
θ ( x, y ) =
π
∞
∑
n=1
(−1)
n +1
n
+1
⎛ n π x ⎞ sinh( n π y Lx )
sin⎜
.
⎟
⎝ Lx ⎠ sinh( n π Ly Lx )
(4)
which is an infinite series solution.
Laboratory:€
For this laboratory, you will write a MATLAB m-file to calculate the temperature from the
analytical solution above at a matrix of x and y coordinates and make a contour plot of the
results. Below are some general instructions on how to write your program. Please refer to the
online MATLAB help files as needed for more complete information. If you have not previously
used MATLAB, or if you have forgotten some of MATLAB’s features or syntax, you should use
the program help material before starting this assignment. In particular, use the “Help” menu
item and “Product Help” sub menu item to bring up the “Help” window. On the left press the
“Contents” tab. Under the “MATLAB” and “Getting Started” submenus read through the
material in the “Quick Start,” “Language Fundamentals,” “Mathematics,” “Graphics,” and
“Programming” sections.
1. Open MATLAB by double clicking on the MATLAB icon on the desktop. Create a new script
(or m-file) using the “File”, “New”, and “Script” menu item. Save the file as “ME554_lab1.m”
(or anything else, but do NOT use a name with spaces) using the “File” and “Save As” menu
item. Recall that a script (or m-file) is simply a list of MATLAB instructions saved in a file with
a “.m” extension. Write a header for your program using “%” for comment lines. Your header
should include the program name, description, your name, date created, date last modified, and a
variable list with definitions (come back and add these later). Continue to use many comment
lines throughout your program to describe each subsequent set of commands and to indicate units
for each parameter.
2. Next, put a line with just
clear, clc
to clear the memory of all variables and clear the screen.
2
3. Set the values for the geometry and boundary condition parameters (Lx, Ly, T1, and T2). You
can do this generally in three ways: simple assignment statement, by reading input from the
screen during execution, or by reading data from a file. Use the first method for some of the
input as follows:
T1 = 100; % deg. Celsius
where the semicolon suppresses output during execution and the units are included as a
comment. Also, try using the second method with an input statement:
T2 = input(‘Input temperature at y = Ly in deg. C:
‘)
4. Create vectors of x and y values of the coordinates at which you need to evaluate the
temperature. First, set Nx and Ny to the number of increments in each direction. Again, you can
do this with any of the methods given is Step 2. Second, calculate the increments between
calculations:
dx = Lx/Nx;
dy = Ly/Ny;
% m, increment in x-direction
% m, increment in y-direction
Third, create the vectors of x and y values using the following:
x = [0:dx:Lx];
y = [0:dy:Ly];
If you do not know what these lines create, enter them in the “Command Window” in MATLAB
with numbers substituted for the variables to test their output. Note that these lines demonstrate
some of the power of MATLAB in that you can create a vector of data in a single line, unlike
FORTRAN or c++.
5. Initialize a matrix to store the data for θ during calculations using the zeros function and set
θ = 1 at your top boundary using:
€
theta = zeros(Nx + 1, Ny + 1);
theta(:,Ny+1) = 1;
This is required for two reasons: (1) to allocate the correct amount of space in memory for your
array and (2) to initialize θ to zero for the summation in the next step.
6. Calculate the value of θ at each (x, y) location. To do this, use two nested for loops to cycle
through each (x, y) location with the following syntax:
for j = 2:Ny
for i = 2:Nx
% insert Equation (4) here
end
end
3
You will need to use a third nested for loop to sum up the terms in Equation (4). Use “…” to
continue your equation on the next line to make it easier to read. Note that for MATLAB “pi”
will give you π with sufficient precision. For this sum, let n go from 1 to 199. Because you are
stopping the calculation at a finite number of terms your results will have truncation error in
addition to machine precision error for double precision calculations. Note that in Equation (4),
for even values of n, the summation term goes to 0 so you can skip these terms to speed up the
calculation. One way to do this in MATLAB is to use “for n = 1:2:199” which increments
n by 2 for each loop instead of the default increment of 1.
7. Convert the calculated θ matrix to temperatures using the variable name “T”.
8. Print out temperature data to the screen for review using fprintf to produce formatted
output. The following lines work well:
for j = Ny+1:-1:1
fprintf('%7.1f', T(:,j))
fprintf('\n')
end
These lines print out the temperatures in the order you would expect for a Cartesian coordinate
system (rotated counterclockwise 90˚ as shown in Figure 2). The ‘%7.1f’ sets the format for the
temperatures as floating point with seven spaces and one number after the decimal place. The
‘\n’ adds a carriage return between lines. Alternatively, you can use (fliplr(T))’ or
rot90(T) to rotate the T matrix before printing it out.
9. Make a contour plot of the resulting temperatures using the following lines:
dT = (T2 – T1)/Nc;
v = T1:dT:T2;
colormap(jet)
contourf(x, y, T', v)
colorbar
Tmax = max(T1, T2)
Tmin = min(T1, T2)
caxis ([Tmin, Tmax])
axis equal tight
title('Contour Plot of Temperature in deg. C')
xlabel('x (m)')
ylabel('y (m)')
where Nc is the number of contours (whose value needs to be set), dT is the step in temperature
between contours, and v is a vector that sets the temperature levels for the contours. The
MATLAB built in function colormap sets the colors used for contours, contourf creates the
filled contour plot where ' after T transposes the matrix, colorbar adds a scale to the plot
where caxis sets the limits, axis equal makes the x and y axes have the same length, and
xlabel and ylabel add x-axis and y-axis labels to the plot.
4
10. Try changing the variables or options for Steps 8 and 9 to alter the output. You can also
explore some of the many other options for plotting data. For example, changing the 4th line to
contourf(x, y, T', v, 'LineStyle', 'none')
will remove the lines between the contours.
Assignment
Submit your lab report as a single pdf file using PolyLearn that contains the items listed below.
1. Set the following values in your code: Lx = Ly = 1 m, T1 = 0 ˚C, and T2 = 100˚C. Include tables
for temperature values at each location for Nx = Ny = 5 and Nx = Ny = 10.
• Label them “Table 1” and “Table 2” along with a descriptive caption above the table.
• Include x and y locations for each temperature on your table so you can determine which
temperatures are evaluated at the same location.
• Use 4 decimal places for temperature data evaluated using double precision.
2. Verify that your solution has converged to at least four decimal places for temperatures on the
interior of your domain. Discuss how limits in machine precision and truncating the infinite
series solution to a finite number of terms influence the accuracy (be quantitative) of your
calculated temperatures.
3. Compare the values in Table 1 and 2 for each case at the same coordinate locations. Are the
temperatures the same or different (be quantitative) at each location where they overlap? Briefly
explain why this result makes sense.
4. Include contour plots for the two cases in step 2. For each figure include “Figure 1” and
“Figure 2” along with a descriptive caption below the figure. Make sure that if you plot or print
out your results in grayscale you use “grayscale” for your colormap so that each temperature has
a distinct shade.
5. Do the contour plots change for these two cases? Briefly explain why this result makes sense
and how this corresponds to the accuracy of your solution.
6. Include a copy of your final m-file.
5
Cartesian Coordinates:
Ti, j = T(x i , y j )
y
€
where
and
and
Δx = Lx N x
Δy = Ly N y
€
€
∆y
€
x i = (i −1) Δx
y j = ( j −1) Δy
€
T1, Ny+1
TNx+1, Ny+1
Figure is
shown with
Nx = Ny = 3
Ly
∆x
T1,2
T2,2
T1,1
T2,1
TNx+1, 1
x
Lx
Matrix Notation:
⎡ T1,1
T1,2  T1,N y +1 ⎤
⎢
⎥
T2,1 T2,2
 ⎥
⎢
[T ] = ⎢ 

 ⎥
⎢
⎥
⎣TN x +1,1   TN x +1,N y +1⎦
€
j (column number)
NOTE: The T matrix is rotated
clockwise 90˚ relative to the
Cartesian coordinate system.
i (row number)
Figure 2. Schematic diagrams of Cartesian coordinate system and matrix notation.
6