Hover over the menu bar to pick a physics animation.

1D Heat Flow with Time Dependent Boundary Condition

Introduction

For diurnal heating due to the sun, it is very important to understand the effects of heat capacity, thermal insulation, and insulator thickness on temperature. There are two effects:

1. The reduction of peak to peak temperature swings as a result of increasing the thickness of an insulator or increasing the heat capacitu of the media.

2. The time lag between the peak solar heating and the peak temperature inside an insulated 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. I found a solution to differential equation 1 as a Green's Function in Wikipedia. The Green's function for time dependent boundary conditions is:

`G(x,t)=x/sqrt(4pialphat^3)exp(-x^2/(4alphat))\ \ \ \ \(2)`

which is a solution for a zero width spike (delta funtion) at `t=s` but can be extended to the sinusoidal time dependence at `x=0` boundary condition that we have here.

`T(x,t)=int_0^tx/sqrt(4pialpha(t-s)^3)exp(-x^2/(4alpha(t-s)))h(s)ds\ \ \ \ \ (3)`

where `h(t)` will be the time dependence at `x=0`. A numerical integral over s appears at first to have divide by zero problems but, as long as `xgt0` and `sltt` by a very small value, the exponential in the numerator goes to zero much faster than the denominator so the integrand goes to zero before `s->t`.

To apply the boundary condition at x=0, we set

`T(0,t)=T_Bcosomega t`
so that:
`h(s)=T_Bcosomega s\ \ \ \ \(4)`

where `T_B` is the amplitude of the temperature swing at the `x=0` boundary and `omega` is the radian frequency of the driving temperature. Or we can just make the boundary temperature constant at `T_B` and observe how the temperature gradually becomes `T_B` across the entire domain. This program uses simple numerical integration of equation 3 with `h(s)=T(0,s)`.

Animation Controls and Graphics

I provide adjustable Diffusivity, `alpha`, driving frequency, `omega`, and distance of thermometer, `x` from the left bound. I also provide radio buttons to choose between sinusoidal driving temperature or constant driving temperature. For finer control of the rate of time increase, I provide a Single Step Checkbox. The thermometer is a simple mercury thermometer that moves up and down with sine input or just up with constant boundary conditions. A checkbox provides a choice to make a color coded temperature plot.

Time Lag of the Temperature

As mentioned in the Introduction, larger distances from `x=0` result in reduced peak to peak temperature range. Also mentioned in the Introduction, was the time lag between the driving temperature and the maximum and minima seen at the thermometer. Heat diffusion is a random walk process so we expect that the lag will increase as the square of the distance, `d_(th)`, of the thermometer. An approximate equation for this lag is

`t_(lag)~=d_(th)^2/(4 alpha)\ \ \ \ \(5)`

where equation 5 assumes that the exponential argument in equation 2 overwhelms the linear `x` and the `t^(-3/2)` factors so the lag time is when the exponential equals `1/e`. Since equation 4 shows that at `t=0` the value of the driving boundary conditon is a maximum, equation 5 refers to the arrival time of the first maximum. Subsequent maxima will arrive at the therrmometer at time `tau=(2pi)/omega` later, where `tau` is the period of the driving temperature.

Solutions for the Envelope of the Extrema

First we have to integrate 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 maxima. 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))\ \ \ \ \(7)`

If we define the scaling parameter `lambda`, equation 7 looks simpler:

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

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

The equation 8 expression for `lambda` is a very important result because it contains the parameters for the reduction of temperature variation with distance from the driving temperature boundary.

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.