# Neumann Boundary Conditions for Heat Equation

I wrote a program that uses the FEM to approximate a solution to the heat conduction equation. I was lazy and wanted to test it, so I only allowed Neumann boundary conditions (I will program in the Dirichlet conditions and the source terms later).

When I input low values for the heat flux, I seem to get reasonable results, as you can see in this http://i222.photobucket.com/albums/dd24/RabidUndeadMonkeys/Neumann_Low.png" [Broken].

When I increase the flux by just a little bit (about 1%), though, it all goes http://i222.photobucket.com/albums/dd24/RabidUndeadMonkeys/Neumann_Chaos.png" [Broken]

In these pictures, red is hot (obviously!), and maximum red value is normalized so that the hottest temperature is pure red. What you can't see is that the temperatures in the first image were very close to zero. The temperature in the second image is about a constant 2 billion. I checked to make sure that the result I got was indeed a solution to Kx = F, (since I am using the conjugate gradient method, it would fail for a non-positive definite matrix), and it was.

My question is this: Is this result actually supposed to happen? Is it that the problem is ill-posed for certain Neumann conditions? Or if it is not, what is likely the problem (I know it's a long shot, but maybe someone has seen this problem before)?

I am using 3x3 Gaussian quadrature for the values of K, and a 2 point Gaussian Quadrature over each side in the elements for the values of F. The values of K(i,j) and F(j) were calculated using the functions given in "The Finite Element Method" by Thomas J.R. Hughes.

UPDATE::
I looked at my stiffness matrix and realized it has exactly one zero eigenvalue, so it's not positive-definite. That is a problem!

Last edited by a moderator: