5

I am attempting to calculate temperature of section of rock in the earth as a function of vertical position in the rock and time. Along with it I am calculating the heat flow through the rock as a function of vertical position in the rock and time. Note that this is a one-dimensional examination of the problem- imagining that the rock is uniform in every direction but up-down.

Please be patient with me: I currently do not have the vocabulary nor the background to talk at the level of familiarity with FEM and differential equations I see in other questions on this site. I am a software developer with some physics background, which has allowed me to get this far.

Description of methodology: I am currently splitting up the section of rock into horizontal slices (depth steps), and then calculating the temperature and heat flow of each depth step at specific time-steps. I am assuming that heat only flows upward from below the section of rock.

As I understand it, the temperature difference between the top and bottom of a depth step would be expressed as: $$ \Delta T = \frac{Q\Delta Z}{\lambda} $$ Where $Q$ is the heat flow at the middle of the depth step in $mW/m^2$, $\Delta Z$ is the thickness of a depth step at this time, and $\lambda$ is the thermal conductivity of the rock (constant).

Additionally, the change in the heat flow between the bottom and top of a depth step could be expressed as: $$ \Delta Q = - C_p \rho \Delta Z \frac{\Delta T}{\Delta t} $$

where $C_p$ is the heat capacity of the rock (constant), $\rho$ is the density (constant), and $\frac{\Delta T}{\Delta t}$ would be the change in temperature of the depth step over the course of this time-step.

As a starting point for the calculations, I know the temperature at the top of the rock layer as a whole (but not the temperature at the bottom of the rock layer). I also know the heat flow value just below the rock layer (but not the heat flow value at the top).

The key piece of complexity to this is that the rock layer is increasing in thickness through time (aka deposition), so that I am adding to the $\Delta Z$ of the top-most depth step and then adding more depth steps to the top of the section of rock as time progresses. As a consequence, individual depth steps get more and more buried, and thus hotter through time, which thus causes the heat flow to decrease as you move upwards in the rock layer (as per equation 2 above).

My method of calculating the temperatures and heat flows for each time-step is as follows:

  1. Increase the thickness of the rock layer, adding depth steps to the top as necessary.
  2. Calculate the temperatures at the top and bottom of each depth step going down the layer (since I know the surface temperature, and thus the starting point), using the heat flow values from the previous time-step.
  3. Calculate the heat flow values at the bottom and top of each depth step going up the layer (since I only know the heat flow at the very bottom of the rock layer), using the temperatures of each depth step calculated in step #2, along with the temperature in the previous time-step.
  4. Repeat steps 2-3 until the solution converges. Note that in step #2 I can now use the heat flows calculated in step #3.

The Problem: This solution does not always converge. As the rock layer gets thicker and thicker through time, each time-step takes longer and longer to converge, until it reaches a tipping point and starts to diverge. What can I do to force this to converge?

I am suspicious that the reason for divergence is because I'm forced to calculate temperature down the depth steps and heat flows up the depth steps. This means that the bottom-most depth step's temperature is sensitive to ALL the changes in temperature above it as per the first equation above. But this then causes the heat flow at the bottom to change from the previous iteration by a fair amount, which then propagates all the way to the top of the rock layer via the second equation above, which then affects the temperature at the top strongly, etc.

Warning: The above problem is actually a major simplification of a larger problem I'm trying to solve- the constants in the equation above are actually temperature dependent and depth dependent in a complicated, non-continuous way. However, the problem I am facing still exists in the simplification as described above.

skybluecodeflier
  • 151
  • 2
  • 1
    The problem we have here is that we can't tell if there's a problem with your methodology (which appears to be fine, at least in the verbal description above), or if there's a subtle bug in your implementation (more likely), or if that's a specific "feature" of your system. It could be an issue with your boundary conditions, since you should have, in principle, a condition that at "infinite depth," your temperature should not exceed some value (the highest temperature in the Earth's core, perhaps)? If it does, your model has a defect that must be fixed. – aeismail Jun 09 '13 at 15:03
  • I would set up a very trivial test case (think 3 elements + BC's), and take some small number of time steps making sure things move exactly as you expect. What is the temperature of the rock layer you add?

    Steady state heat flow is in fact elliptic as you have discovered (all values depend on all boundary conditions). To see this formulate a Laplace/Poisson equation on Temperature.

     – meawoppl Jun 09 '13 at 18:12
  • @aeismail: I have gone over the implementation with a fairly fine toothed-comb, and it seems to give correct results when it is converging. My solution is not exceeding any sort of max temperature at the deepest depth until well after divergence begins. How would I enforce such a limit anyways, given my algorithm? Also, concerning the other comments, I guess my intention here was to rule out the 1st and 3rd possibilities- problems with methodology and divergance being a "feature" of the system. – skybluecodeflier Jun 10 '13 at 19:47
  • @meawoppl: I have done so, and everything appears to be calculating correctly. Note that in my case, the heat flow is not steady state for the simple reason that I am adding material at the top. This changes the temperature of the depth steps through time. This makes the whole system non-elliptic (if I understand the definition of elliptic correctly from Wikipedia). Note the delta-T in the second equation above is the change in temperature through time, not through depth. – skybluecodeflier Jun 10 '13 at 19:54
  • @skybluecodeflier: The boundary conditions would in principle be an "error check"; the algorithm should fail if the condition is violated. – aeismail Jun 12 '13 at 15:53
  • @skybluecodeflier: And another thing to think about—this may be a discretization issue. That is, the way in which you add these "depth steps" may be causing your instability. The more sudden the change, the more likely the system is to "blow up" on you. – aeismail Jun 12 '13 at 16:00

0 Answers0