The aim of this assignment is : 

Given a two-dimensional heat transfer problem, and it's solution; to draw the isothermal contour plot using Gnuplot . 

Problem Statement 

A rectangular plate of dimensions a and b is insulated on two sides. One side is in contact with a mixture of ice and water and hence is at a temperature of 0deg.C. The side BC is subjected to a sinusoidal temperature condition. (refer figure). The value of a is 2m and that of b is 1m. 

(com11.gif) 

The governing differential equation is Laplace's equation, 

(Computer Assignment 12.gif) 

This can be derived from the heat diffusion equation in 2 dimensions for the conditions of steady-state conduction and with no 
heat generation. 

The analytical solution for this problem is given by : 

(Computer Assignment 13.gif) 

For a rectangular plate of length a=2m and breadth b=1m the temperature distribution is given by: 

(Computer Assignment 14.gif) 

Assignment 

1. Divide the rectangular domain (the plate) into a grid. Divide the length into 10 equal parts and the breadth into 5 equal parts. Calculate the temperatures at each point ,or node, of the grid. 

2. Using the data generated above, plot the isotherms, i.e., the lines of constant temperature. You should use Gnuplot for generating the contour plot. 

3. What is the slope of the isotherms at the insulated boundaries? Account for the slope. Hint : There is no heat conducted across the insulated boundary. 

Using Gnuplot to Draw Contour plots 

MEEG 302 

This menu is offered on the basis that you have an account on and are somewhat familiar with the UNIX system on Spencer or Strauss. 

1. Drawing environment : You should have access to a UNIX system such as Spencer (me.udel.edu) or Strauss (strauss.udel.edu) or Brahms (brahms.udel.edu). 

2. Preparing the data file : You should prepare the data file in the block form as shown : 

x     y      T 

0.2  0.2 
0.4  0.2 
0.6  0.2 
0.8  0.2 
............. ............... 
2.0  0.2 
0.2  0.4 
0.4  0.4 
0.6  0.4 
............. ............... 
0.2 0.6 
0.4 0.6 
0.6 0.6 
............. ................ 

Notice that the y values are constant in each block of data while the x values change till the length is covered. The temperatures are calculated for each (x,y) pair using the expression for temperature distribution which has been given. 

You can write a FORTRAN program to calculate the temperature distribution and store it in the data file in the desired block form. Your code can be of the form as follows : 

c Output data for gnuplot 
open (30, file='temp.dat',status='unknown) 
do 200 j= 1, noofydivs 
do 100 i= 1 , noofxdivs 
c your code to calculate the position and temperature x(i),y(i),z(i,j) 
........ 
........ 

write (30,50, x(i),y(j),z(i,j)) 
50 format(3f10.4) 
100 continue 
c Block the data for gnuplot 
write(30,'()) 
200 continue 
.................... 

3. Using Gnuplot to draw the contour plot. 

Here is the plotting code you can use : 

set term vttek # sets terminal type 
set para # changes the meaning of 'splot' from normal #functions to parametric functions. 
set contour # enables contour plotting. 
set nosurface # Whenever 'nosurface' is set, no mesh will be #drawn 
set view 0,0,1 # sets point of view of the plot. This controls the #rotation angles along a virtual 3D coordinate #system aligned with the screen such that the #screen horizontal axis is x, vertical axis is y and #z is perpendicular to the screen. 
set cntrparam bspline 
set cntrpara points 10 
set cntrpara order 10 
set cntrpara levels 15 #sets the number of contours, here it is 15. 
set clip one 
set key -0.3,2.6 #set the contour legend at a specific location 
set xra [0:2.0] #sets the X axis range 
set yra [0:1.0] # sets the Y axis range 
set xla 'X (m)' #sets X axis label 
set yla 'Y (m)' #sets Y axis label 
splot 'num3.f' w l #Plots 3-D surface plot using data from the file #'num3.f' and shows on the screen 
pause -1 #Waits until you press 'enter' before proceeding 
set term postscript "Helvetica" 18 #sets the terminal type to Postscript with #the font set to Helvetica 18pt 
set output "num3.ps" # sets the output Postscript file 
replot #plots again and stores the plot as a Postscript #file,num3.ps 
exit 

Notes: 

1. You can write the commands in a file name with subscript '.gnu' and then type load '<the command file name>' to input this command file. To use Gnuplot type 'gnuplot' at the prompt. Otherwise type 'gnuplot <command file name > ' at the prompt. 

2. The above program assumes that your data file is named as num3.f. You can change it to any other name such as 'temp.dat' or 'num.dat' or 'data.txt' that you might give to your data file. You must use the same name in the command program. 

3. About terminal types : 'vttek' is the terminal type if you are using your account from a Macintosh. If you are using X Terminals then the command should be 'set term X11'. 

4. There are two methods to output your plots, both of which are used in the above program : one is output to a screen which is why you need to set the terminal type; the other is output to a Postscript (.ps) file which can then be printed to a printer later. 

Numerical Solution of Two-Dimensional, Steady-State Temperature Distribution in a Rectangular Domain 

Purpose 

The purpose of this assignment is to apply the finite difference numerical technique to solve the steady-state heat conduction equation, given by 

(comp.proj.21.gif)  (1) 

where T is the temperature, x and y are the physical coordinates of the domain, g is the volumetric rate of internal heat generation (in W/m3 or BTU/ft3*hr), and k is the thermal conductivity. Note that such a heat generation would, for example, be caused by passing electric current through the material. Also note that a common mistake is not to include the heat generation rate in the proper units (recall that it represents energy per unit volume ). 

Upon completion of the general programming, The goal will be to perform three types of analyses: 

* Choose progressively smaller values of [[Delta]]x and [[Delta]]y and observe the behavior of the solution. If the problem has been correctly formulated and solved, the nodal temperatures should converge as [[Delta]]x and [[Delta]]y become smaller. This is one way to check the accuracy of the solution. Analytical solutions are of limited utility in checking the accuracy of a numerical model because most problems which will need to be solved by numerical methods usually do not have an analytical solution. 

* Remember, a converged numerical solution does not ensure an accurate solution to the physical problem. In many cases, the final solution is in serious error simply because the problem was not formulated correctly at the start. (For example, the boundary conditions may be incorrectly imposed). No computer or convergence criterion can correct this type of error. One way to check for formulation errors is to perform some sort of energy balance using the final solution. Hence, an energy balance should be performed as a check on problem formulation. 

* After you are satisfied that the numerical model adequately represents the physical situation, conduct numerical experiments by varying the heat generation rate and thermal conductivity in order to study how they effect the resulting temperature field. 

Graphical representation of the isotherms (lines of constant temperature) may be produced using an already available program called Gnuplot. This program is available on the departmental Unix computer, and was described in the first computer project. 

Problem Definition 

A two-dimensional solid (the third dimension is insulated or is so long that the temperature does not vary in that direction) generates heat internally. One surfaces of the solid is insulated, one is contact with a body at constant temperature, one has constant heat flux input and one surface is cooled with a coolant. We are interested in the steady state temperature distribution. 

As you know, the solution of Equation (1) gives the distribution of temperature in a cartesian domain. This expression is a second-order, linear partial differential equation. As in analytical solutions, we need to specify boundary conditions in order to solve the problem. In this assignment, we will use the mixed type of boundary conditions. 

Let us consider the domain shown below in Figure 1. The block is assumed to have a height H and a width W. A two-dimensional analysis assumes that the domain is either infinitely deep (such as a long bar, with temperature independent of the z direction), or that it is a finite block where the two surfaces parallel to the x-y plane are perfectly insulated. 

(comp.proj.22.gif) 
Figure 1. Problem definition. 

The boundary conditions are: 

T=T1 at x=W 
f([[partialdiff]]T ,[[partialdiff]]x) =0 at x=0 
f([[partialdiff]]T ,[[partialdiff]]y) = qo at y=0 
q= h ( T -Tamb) at y=H (2) 

Numerical Solution Parameters 

Divide the domain into m-1 parts in the x-direction and n-1 parts in the y-direction. Discretization of the domain requires m nodes by n nodes for the mesh size, requiring a dimension statement for the temperature of real T(m,n). The individual nodes are numbered (i,j), where i=1,2,3,.....m and j=1,2,3,....n, related as shown below in Figure 2. 

(comp.proj.23.gif) 
Figure 2. Grid cell notation. 

Note that the steady-state temperature distribution is independent of the thermal conductivity in the absence of internal volumetric heat generation. However, we typically need to assign a value for k to avoid a division by zero error. Such details are important in numerical analyses. For the present assignment, consider the following cases: 

case 1: m=10, n=10, k=10 f(W,m C), g=90 KW/m3 
case 2: m=40, n=40, k=10 f(W,m C), g=90 KW/m3 
case 3: m=10, n=10, k= 100 f(W,m C), g=0 
case 4: m=40, n=40, k=100 f(W,m C), g=0 

In addition, we will use the following for the width, height, and boundary temperatures 

W = 30 cm 
H = 30 cm 
T1 = 100 oC 
Tamb = 20 oC 
h = 500 f(W,m2 C) 
qo = 5000 f(W,m2) 

The high thermal conductivity corresponds to materials such as Brass, while the low value corresponds to materials such as Nickel. 

SUBMIT 

(1a) Equations used for boundaries, interial nodes and corners 
Using energy balance, state the equations used to find the temperatures at (i) corners (ii) boundaries (iii) interior. 

(1b) Program Listing 
A title page, indicating the student's name and the project number, as well as a single listing of the program must be submitted. 

(2) Program Output 
For each case, a graphical output of the isotherms is required. This is achieved by using the program Gnuplot, as explained in the Gnuplot handout. 

(3) Discussions 
Discussions comparing the results of case 1 to 2, and case 3 to 4, should be included. The discussion should focus on the convergence of the numerical solution and the effect of changing properties on the temperature distribution. Will any of the materials start melting ? The total discussion should be confined to one page. 

(4) Formulation Error 
The relative error in the numerical solution can be found by considering a heat balance based on the formula 

(comp.proj.24.gif) 

By definition of steady-state, this should be zero for an exact solution. Note that this is only a relative error since neither the heat input to the block nor the heat out from the block can be calculated with certainty! Calculate this error for all four cases. 

Notes and Hints 

When constructing the program, you may assume that [[Delta]]x=[[Delta]]y. By substituting for the derivatives in Equation (1) with central finite difference expressions and rearranging, we can explicitly solve for the temperature at any node (i,j) as 

(comp.proj.25.gif) (6) 

Note that the temperature at each node is expressed in terms of the temperatures of the neighboring points. To find the temperatures, we need to calculate them at each node several times by sweeping the calculation domain during each iteration. Thus, the solutions improve with each iteration. To start the computations, each node is assigned a temperature (if you don't do this, the computer will set everything to zero, which requires more iterations than a more educated method). One suggestion is to use an average of the boundary temperatures. This is not an initial condition! It is only an initial guess. The difference between an initial guess and an initial condition will become more clear in a later assignment. Also note that each iteration sweep corresponds to an effort to improve the values; it does not represent a step on a time scale. 

The iteration process is stopped when all nodes pass a convergence test. To do this, one must compare values from the current iteration with values from the previous iteration. When the change in value between iterations for each and every point is less than some prescribed value or percentage, the iteration loop is exited and the current values represent the numerical solution. 

In solving for the temperature distribution, one does not need to know the actual physical locations of the nodal points, only the distances between them. Gnuplot requires this information in the form of a specially formatted file (see Gnuplot handout). Fortunately, this information is easily computed with a pair of nested DO loops as follows: 

deltax = W / real(m-1) 
deltay = H / real(n-1) 
do 2 j=1,n 
do 1 i=1,m 
X(i,j) = real(i-1) * deltax 
Y(i,j) = real(j-1) * deltay 
1 continue 
2 continue 

These coordinates can then be written to the mesh file as required by Gnuplot. One final note: as a word of advice, it should be pointed out that in constructing the computer program, one should take care to make it as general and flexible as possible. Therefore, values such as the boundary temperatures are best implemented as variables read through input statements. A small part of the project grades will be based on quality and flexibility of the programs. 

HELP! 

In addition to office-hours, special consulting will be provided by Roopesh Mathur in SPL 227 (March 30-April 4) Sunday-Thursday (7 to 10 pm) in addition to the office hours. You may also use the e.mail facility to send your questions to me or the TA