Finite Element, Time-Dependent heat conduction

[dinaharchery]

Hello,

I am trying to build a C++ program that uses the Finite Element Method to compute the heat conduction/transfer of a set of 3-node triangular elements (2-D) at different time. I am using the explicit forward difference (Euler) and I believe I have everything properly put together except for two sections - initial node temperatures, and boundary corrections at the end of each time-step.

Can anyone explain to me the proper way to initialize the nodes, and how to "correct" boundary nodes at the end of each time step? Currently I initialize all nodes to zero and the boundary "corrections" are done by setting the heat-rate vector to the predefined set of boundary temperatures (arbitrarily set to 2.0).

The algorithm I am using follows:

(1) Form LOCAL 'c' (capacity matrix) and 'k' (stiffness)
(2) Initialize nodal temperatures at time 0
(3) Select time step
(4) Assemble GLOBAL 'C' matrix and GLOBAL 'K' (stiffness)

- At EACH Time Step -
(1) Form ELEMENT heat-rate vector 'r'
(2) Assemble EFFECTIVE heat-rate vector 'R'
(3) Solve for nodal temperatures
(4) Correct nodal temperatures for boundary conditions

I hope someone out there will be kind enough to help me. Thanks in advance.

 

[davefitz]

well, its only been 32 yrs since I took such a course in college.

The theory is worked out in "analytical methods in conduction heat transfer" by glen myers, and the numerical methods ( focusing on minimizing storage space, a problem in the 1970's no longer relevant today) are outlined in "the finite element method in partial differential equations" by Mitchell + Wait.

Why not simply use a freeware solution?

 

[corus]

I'm not sure an explicit scheme is wise as you'll probably need very small time steps to get a reliable solution.
As for the intitial temperatures at time 0? That's for the user to input. The same applies to boundary coditions where the user either sets a Neumann type or Dirichlet condition, if that's how you spell them

Tata

 

[dinaharchery]

Thank you very much for all the replies.

I would use freeware solution but I am building a program that I eventually hope to port to using a GPU as a coprocessor. I know that the explicit solution uses a lumped mass matrix and as such is computationally low but I will move to implicit solutions after finishing the explicit. I know it seems to be a long way around but I want to learn the technique well.

corus, you are right about the user entering the initial temperatures and the time step needing to be small. Am I correct in assuming that the corrected boundary values are placed directly in the solution vector? For example, if node 1 has computed 7.4 and the boundary value of node 1 has been set at 50 then the value of node 1 is re-set to 50 before the next time step.

Once again, thank you.

 

[corus]

I'm not sure what you're trying to do here. If a boundary node has been fixed at a set temperature then it stays at that temperature throughout. If a surface condition has been prescribed as a boundary codition then the nodal temperature will change. I'm not sure what you mean by 'corrected'?

Tata

 

[IRstuff]

Your example indicates serious problem in your calculation flow. A fixed boundary condition should be used in the calculation for determining the heat flow within an element. You should not be in a position to get a divergence on a fixed boundary value.

Have you looked at open source code? If there is one, that would be better than cranking your own unverified code.

TTFN

FAQ731-376

 

[Dave41A]

I'm with IRStuff. I don't follow your statement that " if node 1 has computed 7.4 and the boundary value of node 1 has been set at 50 then the value of node 1 is re-set to 50 before the next time step"

You should NOT be computing node 1. If the temperature of node 1 is in fact specified, then that temperature is known for all time steps.

In other words, if you are trying to solve the linear system [A][x]=, where lower case indicates a "vector" and [A] is a square matrix, you will find that for some nodes, you know the value of "b" and for some nodes you know the value of "x."

The way this is done in practice is that the system of equations is partitioned. I suggest performing elementary row operations on your linear system of equations until all the equations for which x is known are at the top, and all those for which b is known are at the bottom. The matrix [A] is typically sparse, allowing you to form two smaller systems [A'][x'] = [b'] and [A"][x"]=[b"], where primes and double primes indicate sub-matrices of the appropriate matrices. For the sub-system for which the values of x are known, the values of b are then computed directly. For the sub-system where the values of b are known, you either have to invert the sub-matrix A' (computationally intensive) or perform an iterative routine (there are numerous C++ library routines that will do this for you).

Aside from the fact that this is all in C++, the project strikes me as being very a 1970's/80's-ish/Fortran type problem. There's nothing wrong with that, but I would suggest starting with some of the texts from that era (when it was more common to write your own code). More modern texts on FE blow past this stuff. You might even be able to find some Fortran routines that already do all this (that you could call from C++?).

Good luck,

Dave