Finite Difference Heat Transfer

[kevinthenerd]

In undergraduate mechanical engineering heat transfer, one learns about the finite difference method to solve conduction problems for complex geometries with known boundary conditions. These equations build a very sparse linear system. In one dimension, this matrix is tridiagonal. In two dimensions, this matrix is pentadiagonal. In three dimensions, this matrix is heptadiagonal.

I'm trying to write some code to solve a take-home quiz. It's expected that we should use the finite difference method, but I want to be a smartass and use hundreds or thousands of nodes. (If I'm unsuccessful with this method, my backup plan is pretty much completed.) To do this, I can't settle for an n^3 direct solution method, and I can't settle for an n^2 storage method (where n is the number of nodes). I want to store the matrix in a reduced form and use something like the Thomas Algorithm to solve it, but my numerical methods book only covers the Thomas Algorithm for tridiagonal matrices.

I'm running into two problems, though, regarding building the matrix and solving it. I'm trying to solve a two-dimensional problem, and my left hand side matrix takes the form of a three dimensional array. In addition to the X and Y coordinates of the node, I'm using a third dimension (size 5) to store each of the coefficients of the node and its surrounding nodes in addition to an array b to store the right hand side vector.

If I have a node in the middle of a homogenous substance, the 2-dimensional finite difference equations resemble...

1*top + 1*bottom + 1*left + 1*right - 4*middle = 0

This is equivalent to... (assuming Ax=b and assuming the origin on the bottom left of the problem)

1*x[i,j+1] + 1*x[i,j-1] + 1*x[i-1,j] + 1*x[i+1,j] - 4*x[i,j] = 0

This takes the form...

A[i,j,0]=-4 (middle)
A[i,j,1]=1 (right)
A[i,j,2]=1 (top)
A[i,j,3]=1 (left)
A[i,j,4]=1 (bottom)
b[i,j]=0 (right hand side)

How would I go about building and solving such a setup?

I would really appreciate it.

 

[CycloWizard]

When you say "I'm using a third dimension (size 5) to store each of the coefficients of the node and its surrounding nodes", do you mean the thermal conductivity/heat transfer coefficients? It's been a long time since I did this sort of thing (since now I just use finite elements to solve everything real quickly 😛), but I might be able to help with a little clarification.

 

[kevinthenerd]

Originally posted by: CycloWizard
When you say "I'm using a third dimension (size 5) to store each of the coefficients of the node and its surrounding nodes", do you mean the thermal conductivity/heat transfer coefficients? It's been a long time since I did this sort of thing (since now I just use finite elements to solve everything real quickly 😛), but I might be able to help with a little clarification.

Ok, imagine it in 1-D. In a 1-D case, you'd need three coefficients for each node.


1 - 2 - 3 - 4

Eqn1: T(2)+T(0)-2*T(1)=0
Eqn2: T(1)+T(3)-2*T(2)=0
Eqn3: T(2)+T(4)-2*T(3)=0
Eqn4: T(3)+T(5)-2*T(4)=0

In this case, you'd need two dimensions in the array: one for the index of the node and one for the index of its coefficient. In this example, it would be 4 x 3.

Another example:

1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9

This would require an array 9x3


Ok, now for 2-D....

01 02 03 04
05 06 07 08
09 10 11 12
13 14 15 16

In this case, you'd need an array that's 16x5.... 16 nodes with 5 coefficients each. Why 5 each? In 2D, your equations take this form... (example for node 06)

T(07)+T(02)+T(05)+T(10)-4*T(06)=0

 

[kevinthenerd]

In a node that's not on the interior, these coefficients are based on the conductive heat transfer coefficient (k).

 

[CycloWizard]

I think maybe you've overcomplicated things or I'm oversimplifying. You don't need to store the spatial coefficients as you've described, as the central difference (which I believe is what you're using) is always dependent on spatially adjacent nodes. Thus, node ij is always dependent on nodes (i-1)j, (i+1)j, i(j+1), and i(j-1). I have done this in Excel and Basic with no troubles, even for transient problems. I know that the finite element method utilizes a 2-D array to solve such problems in both heat transfer and elasticity, so I know it can be done. I guess I'm still not sure why you would need to store the coefficients as a third dimension in the array. The node 'knows' its coordinates based on its location in the array, even if you don't tell it where it is.

 

[kevinthenerd]

Originally posted by: CycloWizard
I think maybe you've overcomplicated things or I'm oversimplifying. You don't need to store the spatial coefficients as you've described, as the central difference (which I believe is what you're using) is always dependent on spatially adjacent nodes. Thus, node ij is always dependent on nodes (i-1)j, (i+1)j, i(j+1), and i(j-1). I have done this in Excel and Basic with no troubles, even for transient problems. I know that the finite element method utilizes a 2-D array to solve such problems in both heat transfer and elasticity, so I know it can be done. I guess I'm still not sure why you would need to store the coefficients as a third dimension in the array. The node 'knows' its coordinates based on its location in the array, even if you don't tell it where it is.

The node already knows its position, but it doesn't know its temperature. It has to find out this temperature based on the conditions of the surrounding nodes or boundary conditions. Thus, each node has its own equation, and each node's temperature is a variable, giving a square matrix that is very sparse (with a right hand side vector b). The number of rows and columns of this matrix is the product of the rows and columns of the node grid. In other words, if I have a problem set up with 25x15 nodes, giving me 375 equations with 375 unknowns, the resulting matrix has the dimensions 375x375, giving 140625 array elements (most of which are 0). Each row of this matrix has up to five non-zero coefficients, and these coefficients are always based on the surrounding nodes in the original geometry. I'm trying to simplify the storage of this matrix by defining it as a three-dimensional array. The first two dimensions of the array are the x and y coordinates in the original geometry, and the third dimension holds the cardinal direction from the node in question, and the value at the array coordinate holds the equation coefficient.

 

[MoMeanMugs]

Why not use Gaussian elimination?

 

[kevinthenerd]

Originally posted by: MoMeanMugs
Why not use Gaussian elimination?

Gaussian elimination requires N^3 operations and is susceptible to roundoff error. I've thought about this a lot. I'm trying to REDUCE the matrix to a simplified form to overcome this. Even if I don't reduce it, an iterative solution method would be more efficient here.

Example:
A 25x75 node array (not large by modern engineering standards) would take 3.35 million elements of storage using the inefficient method, but condensed it would use only 9375 elements.

If I'm going to use more than a few thousand nodes, I need to optimize the code a LOT better than that.

 

[CycloWizard]

Originally posted by: kevinthenerd
The node already knows its position, but it doesn't know its temperature. It has to find out this temperature based on the conditions of the surrounding nodes or boundary conditions. Thus, each node has its own equation, and each node's temperature is a variable, giving a square matrix that is very sparse (with a right hand side vector b). The number of rows and columns of this matrix is the product of the rows and columns of the node grid. In other words, if I have a problem set up with 25x15 nodes, giving me 375 equations with 375 unknowns, the resulting matrix has the dimensions 375x375, giving 140625 array elements (most of which are 0). Each row of this matrix has up to five non-zero coefficients, and these coefficients are always based on the surrounding nodes in the original geometry. I'm trying to simplify the storage of this matrix by defining it as a three-dimensional array. The first two dimensions of the array are the x and y coordinates in the original geometry, and the third dimension holds the cardinal direction from the node in question, and the value at the array coordinate holds the equation coefficient.

Well, you could do it this way, but I still think you're overcomplicating things. Any problem can be solved using a 2-D system of matrices, and this method is generally more efficient as well (which is why it is the foundation for the finite element method). I'd have to sit down and really figure out how to go about what you're doing, but I don't really have time this week. Maybe next week. 😛

 

[kevinthenerd]

You see, if I build the whole 2-D matrix, it's NOT more efficient. I end up with a matrix of almost all zeroes. Any solution method of that matrix would suck, and the storage would be very inefficient.

By the way, it's not truly pentadiagonal. Each row has up to five non-zero elements, but they aren't arranged in a nice way. Its easy to know where to look for nonzero elements, though, because it's all based on problem geometry,