Hover over the menu bar to pick a physics animation.

Diffusion Equation with Time Dependent Boundary Conditions Solved by Finite Difference Method

Note that this animation is the finite difference equation solution equivalent of the Heat flow 1D with Time Dependent,,, which used classical solutions of the heat flow equation.

In the early days of physics before Ludwig Boltzmann made it clear that heat is really particle motion, heat was considered to be a fluid whose flow depended on temperature gradients. The diffusion equation is what was derived to quite accurately the flow of heat which is the same as the progress that temperature makes when moving from a warmer region to a cooler region. In this animation the temperature swings are shown as a function of conductor (or insulator) thickness. The differential equation that handles temperature diffusion in a solid is:

` (delT)/(delt)=alpha (del^2T)/(delx^2)\ \ \ \ \ (1)`

where `T` is temperature, `t` is time and the diffusivity is `alpha=K/C`, where `C` is the volumetric heat capacity (`("Joules")/(m^3"Kelvin")`) and `K` is the thermal conductivity (`("Watts")/("m-Kelvin")`). A large conductivity results in rapid heat conduction. The insulation parameter, `I=1/K`, is the inverse of the conductivity so the diffusivity in terms of `I` is `alpha=1/(CI)`. Thus a large insulation parameter results in a slow heat conduction. To solve equation 1 Vs `x` an `t` I use the finite difference method. The second derivative is represented by

`(del^2T[i])/(del x^2)~=(T[i+1]-2T[i]+T[i-1])/(deltax^2)\ \ \ \ \(2)`

so equation 1 becomes:

`T'[i]~=T[i]+alpha(T[i+1]-2T[i]+T[i-1])((deltat)/(deltax^2))\ \ \ \ \ (3)`

where `T'[i]` is the new iteration result.   For stability, the quantity `alpha(deltat)/(deltax^2)` must be less than 1/2. Thus we can just set `(deltat)/(deltax^2)` equal to 1/2 and constrain `alpha<= 1`.

Method for Making the Left Boundary Temperature Time Dependent

Since we want to see the effects of changing the input temperature, I made the temperature at the left (`x=0`) boundary change sinusoidally with time:

`T[0,t]=T_Bcos(omegat)`

where `T_B` is the maximum temperature, `omega'` is the driving freqquency (radians/secomnd) and t is time. By itself that specification of temperature distribution does not work with equation 3. Instead I had to use a modification that has a finite width `sigma`in the x direction:

` T[x,t]=T_Bcos(omegat) exp(-(x^2)/(2sigma^2))\ \ \ \ (4)`

where `sigma` is much less than the range of `x`. Then I used this equation to created a matrix of many rows that had sequential values of 't' separated by `deltat`. In the animation, for each frame, a new row of this matrix is added to T and then eqution 3 is applied to the new sum.

Solutions for the Envelope of the Extrema

First we have to iterate equation 3 to numerically determine the envelope. To do that I defined an array of "listeners" equally distributed along the range of the x axis. Besides their x position, the listener objects have minimum and maximum temperature components. Each listener computes the temperature at each time iteration of the animation. When the listener determined a temperature that was higher than its present maximum the listener maximum became that temperature. When the listener determined a temperature that was less than its present minimum the listener minimum became that temperature. A plot of the maxima and minima becomes the envelope of the extrema. This plot stabilizes after a short period of time. It turns out that an exponential Vs `x` fits the numerically computed envelope. That exponential is:

`T_(extrema)(x)=T_Bexp(-x/sqrt((2 alpha)/omega))\ \ \ \ \(5)`

If we define the scaling length `lambda`, equation 5 looks simpler:

`lambda=sqrt((2 alpha)/omega)=sqrt((2(K/C))/omega))\ \ \ \ \(6)`

`T_(extrema)(x)=T_Bexp(-x/lambda)\ \ \ \ \(7)`

Equations 6 and 7 are very important results because they specify the reduction of temperature variation with distance from the driving temperature boundary which is really the effect of adding insulation.

Conclusions

For diurnal solar heating, these results show the effects of decreasing conductivity of walls of a house. The temperature inside is greatly stablized by increasing wall insulation (reducing `K`). In addition, these results show the effects of the very large thermal inertia of a cave. While the walls of the cave do not have particularly small conductivity, the cave wall thickness heat capacity (`C`) is sufficient to allow the temperature inside to fluctuate very little with solar heating. The time lag for a cave is so large that it often takes months of increased solar heating to significantly increase the temperature inside.

Followon Suggestions

Since equation 7 contains the rolloff of the extrema for any single sinusoid, it also pertains for each of a series of sinusoids. That means that it can pertain to a Fourier series expression for something like a square wave of a triangular wave.