ME 314 Heat and Mass Transfer Laboratory

--

The objective of this experiment is to gain familiarity with numerical methods and use them to solve a two dimensional steady state heat convection problem in an enclosure.

 
INTRODUCTION:

In the previous lab we solved a two dimensional steady state heat conduction problem through a solid using the numerical method. More practical situations, as discussed in the previous lab, involve advection heat transfer, as a result of bulk motion of a fluid, in addition to conduction. The combination of advection and conduction in the presence of a solid surface is commonly referred to convection.

You will see that the addition of advection complicates the analysis much more. In the previous lab where we dealt only with conduction , temperature was the only dependent variable. The addition of advection contributes three more dependent variables, the (x) and (y) velocities and pressure. In addition to the energy equation used in the previous lab, we will have to add three more equations, the x and y momentum equations and the continuity equation.

We will first derive the discretization equation for the temperature distribution using the energy equation while assuming that we know the velocity distributions. We will see that the x and y momentum equations are of similar form as the energy equation and can be discretized in a similar fashion.

The velocities can then be calculated using the discretized x and y momentum equations respectfully, the discretized continuity equation, and a pressure field. You will see that the pressure term in both momentum equations will cause some difficulty. A pressure equation must be constructed from the substitution of the momentum equations in the continuity equation allowing one to solve for the pressure field using guessed velocities. An iterative procedure must be used to continually improve the guessed velocities in the pressure equation until the correct pressure is obtained. The correct pressure used in the momentum equation results in velocities that satisfy the continuity equation. Knowing the pressure field, the velocity fields can be calculated. Once the velocity fields are calculated, the temperature distribution can then be calculated.

For this lab we will solve a two dimensional steady state convection problem using a numerical computer program.

 
THEORY:

Let's start by looking at the thermal energy equation under steady state, two dimensional conditions which is Equation 6.44 in your textbook [2]. In addition we will assume that viscous dissipation can be neglected as well as work done by pressure. The corresponding equation is
 

(1)

 
where k is the thermal conductivity, T is temperature, and S is the volumetric heat production within the solid (which is the same as "q dot" in your textbook), and u and v are the x and y velocities.

If we combine the advection and conduction term and call this the total flux (J) for simplicity, Equation 1 becomes,
 

(2)

 
where,
 

(3)

 

(4)

 
The symbol G is the ratio of "k" and "c_p". From Equations (3) and (4) we see that the total flux is a combination of advection and conduction. For now let's assume that we know the flow field (i.e., the velocity components). How we obtain the flow field will be discussed later.

Consider a point P which is in some plane inside the tank in Figure 1. This figure illustrates the five point grid cluster at point P. The volume for control volume in Figure 1 is simply (Dx)(Dy)(1) if we assume unity in the (z) direction. Integrating over the control volume we get,
 

(5)

 
where J_e, J_w, J_n, and J_s are the integrated total fluxes over the control volume faces; that is

over the interface "e".

 
The continuity equation, which is equation 6.25 in your text [2], is
 

(6)

 
We can integrate over the control volume and obtain,
 

(7)

 
 

where F_e, F_w, F_n, and F_s are the mass flow rates through the faces of the control volume. If ru at point "e" is taken to prevail over the whole interface (e), then
 

(8a)

 
Similarly,
 

(8b)

 

(8c)

 

(8d)

 
If we multiply Equation (7) by T_p and subtract it from Equation (5) we can essentially combine the continuity and energy equation together to form
 

(9)

 
Now we would like to put this equation into a general discretized form similar to heat equation in the last lab. We can do so by relating each term in Equation (9) to the coefficient at the face multiplied by the temperature differences between the two grid points,
 

(10)

 
such that a_E is,
 

(11)

 
Some steps have been omitted because of the complexity but can be found elsewhere [1]. The coefficient "a" is somewhat different from what we saw in the conduction problem. Recall last time we used a piece-wise linear profile between the grid points. Here we have made no assumption about how properties behave between the grid points. Instead we use the general function (A|P|). The function ( A|P_e|) defines the profile scheme between the two grid points at the interface "e" and it is a function of the Peclet number (P) which is the ratio of the relative strengths of the advection and conduction terms. This will be discussed later. The function "max" just takes the larger of the two arguments. The reason this function appears in Equation 11 is to ensure that the coefficient remains a positive value. Positive coefficients are necessary for a converging solution. The explanation is complicated and can be found elsewhere [1].

The final two dimensional discretization equations can finally be written as,
 

(12)

 
You might recognize that this is same general discretization form that you saw in the conduction problem.

The linear piece-wise profile scheme unfortunately is not a good assumption for how properties behave between grid points for the convection problem. The reason is very clear if we look at the exact solution for a one dimensional case between two grid points. The exact solution is found by integrating the differential equation between two grid points and is
 

(13)

 
where P is the Peclet number. The Peclet number is very important. This number is the ratio of relative strengths of the advective and conductive terms in Equations (3) and (4) and is defined as
 

(14)

 
Figure 2 shows graphically the solution of T as a function of the distance x between two grid points for various Peclet numbers. From this figure we can see why we cannot use the piece-wise linear profile. This scheme is only valid for small Peclet numbers, typically |P| <2. As the Peclet number tends toward zero we have pure conduction and a linear relationship exists. However, when the Peclet number is large, advection dominates and the profile is highly non-linear. So it is necessary to develop a profile scheme that will best represent the behavior of the properties between two grid points for a wide range of Peclet numbers.
 

So far we have formulated the procedure for solving for the temperature in the presence of a given flow field. Now we would like to know how to obtain knowledge of the flow field. The procedure for calculating the flow fields is quite complicated and only the methodology will be presented. For the x and y velocities we use the x and y momentum equations which are equations 6.34 and 6.35 in your textbook [2], These equations can be simplified and put in a general discretization form similar to the energy equation (Equation 12). They could then be solved in a similar manner as the energy equation. However these equations contain a pressure term which is unknown.

One way of dealing with the unknown pressure field is to combine both of the momentum equations into the continuity equation to develop a pressure equation for which an initial pressure can be solved. After substituting both momentum equations in the continuity equation you will find that the equation can be manipulated into a general discretized form similar to the energy equation and momentum equations discussed earlier with pressure as the independent variable. The source term in this equation will contain the velocities, which will have to be initially guessed. Using the guessed velocity fields in this equation, a first approximated pressure field p* can be calculated. If we use this pressure in the momentum equations we can solve for new velocities v*, and u*.

These velocities usually must be corrected before using them in the pressure equation to help convergence. This can be done by using a pressure correction equation in combination with a velocity correction equation, which will allow a corrected velocity to be calculated. The corrected velocities are then put back in the pressure equation to get a new pressure and after several iterations we would hope to arrive at a solution for a pressure field that is very close to the actual pressure field (i.e. this pressure field supplied in the momentum equations provides velocity fields which satisfy the continuity equation). This is basically a summary of the widely used SIMPLER algorithm developed by Patankar [1].

One last important point should be discussed. The control volumes for the x and y momentum equations are different than those discussed for the energy equation. Figure 4 shows the x and y momentum control volumes. The control volumes are called staggered control volumes. The control volume for the x momentum is staggered in the x direction only, and the y momentum control volume is staggered in the y direction only. As a result the velocities are stored between the grid points as opposed to at the grid points like temperature and pressure.

There are two basic advantages to doing this. The first advantage is the elimination of the possibility for an actual "wavy" velocity pattern to be mistaken as a uniform velocity field. For example, consider a one-dimension case where the continuity equation states that du/dx=0. If the velocities were not staggered our discretized continuity equation, using a piece-wise linear profile for the midway location of the control volumes, would be
 

(15)

 
As a result, the continuity equation demands that alternate grid points u_E and u_W must be equal. For an actual "wavy" velocity field such as
 

Just like velocity, an actual "wavy" pressure field could also be mistaken as a uniform pressure field if the staggered control volume approach for the momentum equations was not employed. However, because the momentum control volumes are staggered, the pressure difference between two adjacent grid points is the natural driving force for the velocity component located between these grid points, and, hence, two adjacent pressures are used instead of two alternating pressure nodes (see Figure 5).

Another advantage is one of convenience. When calculating the mass flow rates for the energy equation at each face, the velocities do not have to be interpolated between two nodes, the velocities are stored at the faces of the temperature and pressure control volumes and can be used directly as they are (see Figure 5).

Many small details have been omitted for reasons of complexity. Only the basics have been presented. If you are interested in the details please refer to the reference [1].

 
COMPUTER PROGRAM:

The computer program that will be used for this experiment, written in FORTRAN, is me314p2.exe^*. The program utilizes the numerical method to solve the temperature and velocity distributions for a two dimensional cross section of a tank consisting of water at steady state. The tank is similar to Figure 1 and is insulated on the top and bottom. The program uses a grid consisting of 26X36 nodes. The (x) and (y) locations of the nodes are stored in an input file, "grid2.dat", for which the user can specify (dimensions in centimeters). During the execution of the program, the user is asked to input data such as density, viscosity, etc., and boundary conditions for the left and right side of the tank. Temperatures must be in Kelvin. After the program has executed the node locations velocities and the temperatures at each node are stored in the output file, "undump.dat". The information from this file can be used in a plotting routine to plot the temperature distribution and streamlines in the x and y direction. A plotting program^** for Matlab is available that reads the data from "undump.dat" and plots a 3-D mesh or contour plot of the temperature distribution. The program name is "graph.m".
 

 
PROCEDURE:
 

1. Obtain a grid distribution from the lab instructor.

 
2. Store the (x) and (y) location values in the input file.

 
3. Execute the program for various conditions as specified by the instructor.

 
4. Plot the temperature distributions (contour plot) and streamlines using a plotting routine or the Matlab program provided.

 
5. Refer to your fluid dynamics text if your need to refresh your memory of streamlines.

 
6. Discuss briefly the following areas:

1. Based on the plots discuss whether or not your results seem reasonable.
2. Discuss some ways that would help you determine if the program solution is realistic.
3. Comment on the accuracy of the numerical method and how it could be improved.

 
REFERENCES:

[1] Patankar, Suhas V., Numerical Heat Transfer and Fluid Flow. Hemisphere, McGraw-Hill, New York, 1980.

[2] Incropera, F. P. and D. P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley and Sons, Toronto, 1990.