ME 540

Numerical Solutions Applied To Heat Transfer And Fluid Mechanics Problems


Approximation of PDEs for Computer Solution

Four basic approaches

1) Control volume methods (CVM)

Complete region is subdivided into control volumes. The nodes are located at the center of the control volumes. Word statement of the conservation equation is used to form difference equation, or the differential form of the conservation equation is integrated over the control volume to form difference equation.

2) Differential equation methods (DEM)

Construct grid network with nodes located on boundaries. Form difference equation using Taylor Series or the exact solution to a one dimensional form of the conservation equation.

3) Spectral Methods

Express all or part of the solution as a Fourier Series expansion. Substitute into equations and solve for coefficients. Troubles arise quickly in terms that are non-linear in the solution variables. We won't cover this approach in this class.

4) Finite Element Method

Express the solution as a parametric functional form. Use expressions minimizing the difference between this parametric function and the actual solution to determine values of the free parameters. We won't discuss this method in this class.
 

Control Volume Method

This method has high value when rigorous conservation of mass and/or energy is very important in the computational simulation of the physical system. This is often the case in transient simulations of internal flows, or thermally isolated systems.

Word statement of conservation equation for a control volume:

Rate of accumulation of a field variable in a volume =

Net rate of flux of the field variable across the surface bounding the volume

+ the net rate of generation of the field variable within the volume
 

The flux of  across the boundary per unit surface area of the boundary

is represented by J. It is composed of an advective and diffusive

component. For energy in a flow in the "x" direction

Example: Steady state conduction in 2 dimensions with internal heat generation

Word statement


 

For a simple conduction problem is no advection of energy.

Thus, for a rectangular control volume centered at (x,y) with height:

Assume flux is uniform over the control volume face. The expressions for flux terms are:


 

Adapting an integer index notation for the volume centers and a half integer index notation for the edges, the final expression becomes
 


 

If you think about application of the above equation at all volumes in the solution region, you will realize that introduction of temperatures at the surfaces results in far more unknowns than available equations. We deal with this through "flux matching". Flux into a surface from one volume must match flux out of the surface into the adjacent volume. For example at the face indexed as i+1/2, j the following equality holds:

or

which gives

For uniform conductivity this gives what you would expect

If one conductivity is much larger than the other what is the answer? What's your physical explanation for that result?

Flux matching is also important for constructing the equations of volumes next to boundaries. Assume surface i+1/2, j has fluid on the other side.

This gives a surface temperature :

Normally we eliminate temperature at volume surfaces using the above internal and boundary equations before solving the coupled set of volume centered energy equations for the volume centered temperature. For n volumes, this gives us a set of n linear equations and unknowns. The trick is to define a 1-dimensional array of unknowns that maps to the spatial 2-dimensional array of temperatures. If we have nx volumes in the x direction, one way to do this mapping is to let k=nx(j-1)+i then define variables in our linear system as:

xk = Ti,j

The linear system has the form Ax=b, where the kth row in A and bk are determined by the energy equation for the corresponding volume i,j.
 

Integration of differential equation over the control volume


 


 
 
 
 


 
 

Assume heat source is constant over the volume and derivatives constant over the bounding surfaces, then we get:


 
 

If you examine the definitions of energy fluxes, you will see that we are now at the same equation obtained from the "word problem" expressing basic conservation in a volume.
 

Implementation of the Solution in a Program
 

Take a look at my Fortran 90 notes for some very basic guidance on creation of a well structured program. To make life easier, I expect you to adhere to good programming practices when you program for a homework. The main program should just be an outline. Create separate subprograms for each important step in the solution (e.g. initialize terms in your linear equations, solve the equations, output the solution). Parameters such as conductivity, volumetric heat source, number of mesh cells, and mesh cell sizes should all be stored in variables that are initialized once in your program. Include adequate comments describing the purpose of each subprogram, defining all subprogram arguments, and describing all key steps in the program. Your goal is to create a program that will be easy to modify as the homework assignments evolve, and will be easy for you to pickup a year or more from now to help solve some new problem.

As a first cut, cast your linear system as a full matrix and use a brute force general matrix solver to get an answer. Take a look at suggestions from my Fortran 90 notes for a linear system solver. LINPACK or LAPACK subprograms are available in both Fortran and C at www.netlib.org. Be sure to test your linear solver before using it in a finite volume calculation. Take a matrix A generated by the finite volume equations, pick a random set of values for vector x, and generate a right hand side b=Ax. Now feed A and b to your solver and see if you get the original values for x back as the solution.

A note on language differences:

My instructions above for assigning elements of xk are good for the following Fortran 90:

k = nx*(j-1) +i
x(k) = T(i,j)
However, in C if i continues to be the x index and j the y index, I would structure my 2-dimensional array as follows:
x(k) = T[j][i]
Why?


Other Lectures / Home


Created by  Frank Schmidt
Maintained by John Mahaffy : jhm@psu.edu