Changing the time step in a 2D transient heat transfer

[rconnor]

Homework Statement

I am working on a Matlab code that will solve a temperature grid for a 2D transient heat transfer problem. The code seemed to be working well but when I changed the timestep (to find timestep independance) I noticed that the larger the timestep, over the same duration, the cooler the temperature field was. Now, a larger timestep should mean that each iteration has a greater change in heat transfer but with less iterations, which should provide, more a less, a similar temperature field (although relatively less accurate than the shorter time step). But I get a much much cooler temperature distribution.

Homework Equations

I am using the implicit finite difference method where I assume the heat flow is into the control volume from the surrounding nodes. The temperature from the node above the node being calculated is called Tn (or Tnorth) the temperature from the node to the right is Te (or Tsouth) etc...The equation (for interior nodes) is:


Tp(new)= Fo*Tn(new) + Fo*Ts(new) + Fo*Te(new) + Fo*Tw(new) + Tp(old) / (1+4Fo)

Where Fo is the Fourier number.


In the formula above, it maybe hard to tell that the time step does not seem to affect the temperature difference but the equation is a rearrangment of an eariler form of the energy equation:

Tp(new)-Tp(old)*rho*Cp*V/dt = k*dy*dz*(Tw-Tp)/dx + k*dy*dz*(Te-Tp)/dx + k*dx*dz*(Tn-Tp)/dy + k*dx*dz*(Ts-Tp)/dy

Here the only place that the timestep (dt) comes up is in the denominator of the Tp(new)-Tp(old) term. I understand way it is there but when I think about it from the perspective of the code, it doesn't seem to effect the value of Tp(new)-Tp(old). Since the k*dy*dz*dT terms are not effected by the timestep, I think this is were the code is falling apart. But thermodynamically I can't see anything wrong with the formula and I know its right because I have used it in school before a lot.

The Attempt at a Solution

What I believe is happening is that nothing in the formula (showen below) dictates that a greater timestep should have a greater temperature difference (due to the increased time which energy has to flow into the node) and thus the temperature distribution is based on the number of iterations in the duration. Since a greater timestep will have less iterations for the same duration, the temperature field is colder.


Is there some reason why this formula cannot be used for the computational analysis or something else I'm missing? I've been going it over and over and over but I can't find anything wrong. I know there is nothing wrong with the logic of the code, it has to do with the timesteps in some way. I believe that the issue is somewhere in the energy balance equation because I can't see where in the formula the length of the timestep effects the temperature difference from the old time to the new time.

If anyone has any insight into this, it would be greatly appreciated! Thank you very much in advance!

 

[Mapes]

I'm not sure where the problem lies, but your equations look fine, and this approach is valid; I've used it myself to simulate heat transfer. I presume you're then solving the equations simultaneously in MatLab to get the temperature field at each time step? Some amount of convergence is to be expected as you reduce the time step, but is it the case that if you half the time step (and double the number of iterations), the temperature change at anyone node is halved? This could indicate that there's an extra $\Delta t$ floating around in the code. Try the usual methods of simplification: reduce the number of nodes to two, reduce the number of durations to one, compare to an analytical solution, etc.

 

[rconnor]

Some amount of convergence is to be expected as you reduce the time step, but is it the case that if you half the time step (and double the number of iterations), the temperature change at anyone node is halved?

Yes, this is exactly what seems to be happening but except reversed, halfing the time step means that the temperature at each node is double (because the output seems to be independant of the timestep, so the only change that halfing the time step causes is doubling the iterations and therefore the grid has twice as many iterations to heat up in). I've been looking through the code and I can't find anything anywhere, but I will continue to search. I am currently retrofitting the code to solve of a very simple grid as you suggested.

The major issue is that the code appears to be independant of the timestep, meaning the only effect it has on the code is altering the number of iterations. The temperature difference between time steps SHOULD increase with decreasing time steps which could balance out the fact that it would have less iterations. But the first part doesn't seem to be working...so any input on why this might be would be greatly appreciated.

Thank you for the suggestion

 

[minger]

How knowledgeable are you with numerical solving?

Firstly, if you're solving a numerical problem, then theoretically, the time step shouldn't matter. In the real world it does, but in a situation like this, you should be able to take a helluva time-step without getting stability problems.

If by changing your time-step, you get a significant change in your answer, then as mentioned before, you have your equations written incorrectly. Sorry. It sucks and debugging is almost always 90% of time spent programming.

If you can read FORTRAN, I have found an old program that I still have around from my Conduction class (assuming this is the class that you're taking too) that you could look at. To give you a snippet, the RHS module interior points look something like: (note: m is the number of points in a "row".

rhsvec(1) = T(1)*(2.0_8 + 2.0_8*(dxdy**2) - 2.0_8*Fo) - T(1+1) - &
            T(1+m)*(2.0_8*(dxdy**2)) - 2.0_8*Tb
DO i=2,m-1
rhsvec(i) = T(i)*(2.0_8 + 2.0_8*(dxdy**2) - 2.0_8*Fo) - T(i+1) - T(i-1) &
            -  T(i+m)*(2.0_8*(dxdy**2))
END DO

And the left hand side looks something like:

! create the interior points,noting that i=1 has a slightly different
! stencil since i-1 is actually at Tb
 matrix(1,1)   = -2.0_8 - 2.0_8*(dxdy**2) - 2.0_8*Fo
 matrix(1,1+1) = 1.0_8
 matrix(1,1+m) = 2.0_8*(dxdy**2)
DO i=2,m-1
 matrix(i,i)   = -2.0_8 - 2.0_8*(dxdy**2) - 2.0_8*Fo
 matrix(i,i+1) = 1.0_8
 matrix(i,i-1) = 1.0_8
 matrix(i,i+m) = 2.0_8*(dxdy**2)
END DO

If you could use some more help (reference), let me know.