**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:

x_{k} = T_{i,j}

The linear system has the form Ax=b, where the kth
row in A and b_{k} 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 x_{k}
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?