Periodic Steady State
=====================

We look at the heat equation again,
but now with periodic inputs.

Systems with Periodic Inputs
----------------------------

A question concerning the decay of periodic inputs is
*How deeply within the earth is felt the temperature at the surface?*

The temperature at the surface is periodic:

* warm (or even hot) during the day,

* cold when the sun sets and at night.

We have a system with a periodic input.
The above question can then be formulated in a more general way:
*How long before the periodic inputs have decayed?*

Some example questions are below:

* Geothermal heat pumps, how deep to bury?

* Heating of concrete, how thick should walls be?

* How deep do worms go to prevent from freezing?

The question we will aim to solve is
*How deep do we have to dig}
to store something at constant temperature?*

The result of a simulation is shown in :numref:`figPeriodicTemperatures`.

.. _figPeriodicTemperatures:

.. figure:: ./figPeriodicTemperatures.png
   :align: center

   A simulation of periodic temperature simulations.

The Helmholtz Phasor Method
---------------------------

Considering the heat equation,
the temperature :math:`u(x,t)`
at depth :math:`x`, at time :math:`t`,
should be a periodic function with decaying amplitude.

Thus we have to solve the heat equation

.. math::

   \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}

with boundary conditions, using idealized numbers:

* :math:`u(0, t) = \sin(\omega t)` is the temperature at the surface,

* :math:`u(\infty, t) = 0`, the surface temperature 
  is not felt at :math:`\infty`.

We covered the method of seperation of variables,
and thus come up with a first proposed solution:

.. math::

   U(x, t) = X(x) \sin(\omega t - \theta(x)),

where

* the amplitude :math:`X(x)` decreases for increasing :math:`x`,

* :math:`\theta(x)` is the phase lag at depth :math:`x`.

The problem with this proposal is that substituting :math:`U(x,t)`
in the heat equation and then solving for :math:`X(x)`
and :math:`\theta(x)` is hard.

Helmholtz proposed

.. math::

   \rho(x, t) = e^{i (\omega t - \theta(x))} X(x)
              = e^{i \omega t} Y(x),

which is called *a phasor solution*.

.. index:: phasor solution

Substituting :math:`U(x, t) = e^{i \omega t} Y(x)` into the heat equation
will lead to ordinary differential equations.

In substituting :math:`U(x, t) = e^{i \omega t} Y(x)` into

.. math::

   \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}

we compute

.. math::

   \begin{array}{rcl}
   {\displaystyle \frac{\partial U}{\partial t}}
   & = & i \omega e^{i \omega t} Y(x) \\
   {\displaystyle \frac{\partial^2 U}{\partial x^2}}
   & = & e^{i \omega t} Y''(x).
   \end{array}

This leads to the ordinary differential equation

.. math::

   Y''(x) = i \omega Y(x).

Let us compute a symbolic solution:

::

   julia> using SymPy
   julia> x = Sym("x")
   julia> omega = Sym("omega")
   julia> Y = sympy.Function("y")
   julia> ode = diff(Y(x), x, 2) - sympy.I*omega*Y(x)
   julia> println(dsolve(ode, Y(x)))
   Eq(y(x), C1*exp(-x*sqrt(I*omega))
          + C2*exp(x*sqrt(I*omega)))

So ``SymPy`` computed

.. math::

   Y(x) = C_1 e^{-x\sqrt{i \omega}} + C_2 e^{+x \sqrt{i \omega}},
   \quad i = \mbox{\tt I} = \sqrt{-1}.

We expect the solution

.. math::

   Y(x) = C_1 e^{-x \sqrt{i \omega}} + C_2 e^{+x \sqrt{i \omega}}

to decay.  Therefore, consider the real and imaginary
parts of :math:`\sqrt{i}`

.. math::

   \sqrt{i} = 
   \frac{\sqrt{2}}{2} \left( \vphantom{\frac{1}{2}} 1 + i \right),

which leads to

.. math::

   Y(x) = C_1 e^{-x 
   \frac{\sqrt{2}}{2} \left( \vphantom{\frac{1}{2}} 1 + i \right)
   \sqrt{\omega}}
   + C_2 e^{+x
   \frac{\sqrt{2}}{2} \left( \vphantom{\frac{1}{2}} 1 + i \right)
   \sqrt{\omega}}.

In our model: :math:`x \geq 0`, and we want decay.
Thus: :math:`C_2 = 0`.

To work with idealized temperatures,
we normalize :math:`C_1 = 1`.

Returning to the phasor solution,
we found :math:`\displaystyle Y(x) = e^{-x 
\frac{\sqrt{2}}{2} \left( \vphantom{\frac{1}{2}} 1 + i \right)
\sqrt{\omega}}` for use in
:math:`\displaystyle  U(x, t) = e^{i \omega t} Y(x)`,
where :math:`Y(x)` appears in the phasor
:math:`\rho(x, t) = e^{i (\omega t - \theta(x))} X(x)
= e^{i \omega t} Y(x)`.

We rewrite :math:`Y(x)` into two parts,
an amplitude and phase lag part:

.. math::

   \begin{array}{rcl}
   Y(x) & = & {\displaystyle e^{-x
   \frac{\sqrt{2}}{2} \left( \vphantom{\frac{1}{2}} 1 + i \right)
   \sqrt{\omega}}} \\
   & = & {\displaystyle e^{-x
   \frac{\sqrt{2}}{2} \sqrt{\omega}}
   e^{-x i \frac{\sqrt{2}}{2} \sqrt{\omega}}} \\
   & = & {\displaystyle 
   \left( \underbrace{e^{-x   \sqrt{\omega/2}}}_{\mbox{decays}} \right)
   \left( \underbrace{e^{-x i \sqrt{\omega/2}}}_{\mbox{oscillates}} \right) } \\
   & = & \underbrace{X(x)}_{\mbox{amplitude}} ~~ \cdot ~~
   \underbrace{e^{-i \theta(x)}}_{\mbox{phase lag}}
   \end{array}

The phasor solution is

.. math::

   \rho(x, t) = e^{i \omega t} Y(x),

with

.. math::

   Y(x) = \left( e^{-x \sqrt{\omega/2}} \right)
          \left( e^{-x i \sqrt{\omega/2}} \right).

Thus, the :index:`phasor solution` is

.. math::

   \rho(x, t) = 
   \left( e^{-x \sqrt{\omega/2}} \right)
   \left( e^{i \left( \omega t - x \sqrt{\omega/2} \right) } \right).


An Application
--------------

Let :math:`\alpha` be the diffusion constant of material.
Introducing :math:`\alpha` in the phasor solution:

.. math::

   \rho(x, t) = 
   \left( e^{-x \sqrt{\omega/(2\alpha)}} \right)
   \left( e^{i \left( \omega t - x \sqrt{\omega/(2\alpha)} \right) }
   \right).

The *depth of penetration* is
:math:`\displaystyle \delta = \sqrt{\frac{2 \alpha}{\omega}}`,
the reciprocal of :math:`\sqrt{\omega/(2\alpha)}`.

.. index:: depth of penetration

If :math:`x = \delta`,
then the amplitude is :math:`e^{-1} \approx 0.3678`.
The depth of penetration is the depth at which amplitude variations
have dropped to :math:`e^{-1} \approx 37\%` of the surface temperature.

Using :math:`\displaystyle \delta = \sqrt{\frac{2 \alpha}{\omega}}`
in :math:`\displaystyle \rho(x, t) = 
\left( e^{-x \sqrt{\omega/(2\alpha)}} \right)
\left( e^{i \left( \omega t - x \sqrt{\omega/(2\alpha)} \right) }
\right)` gives,
if we focus on the real part:

.. math::

   T(x, t) = T_0 + A e^{-x/\delta}
   \cos\left(\omega t - \frac{x}{\delta} \right),

where

* :math:`T_0` is the temperature at the surface, and

* :math:`A` is the maximal variation of the surface temperature.

The parameters for a numerical simulation are

.. math::

   \begin{array}{rcl}
   \alpha & = & 1.0\mbox{E}\!-\!6 \\
   P & = & 24 \cdot 3600 \\
   \omega & = & 2 \pi/P \\
   T_0 & = & 10 \\
   A & = & 10
   \end{array}

For :math:`P`, we take a 24-hour period, measured in seconds.

The evolution of the simulation is shown in :numref:`figTempHour1`,
:numref:`figTempHour4`, :numref:`figTempHour9`, :numref:`figTempHour13`,
:numref:`figTempHour16`, and :numref:`figTempHour21`.

.. _figTempHour1:

.. figure:: ./figTempHour1.png
   :align: center

   The simulation after one hour.

.. _figTempHour4:

.. figure:: ./figTempHour4.png
   :align: center

   The simulation after four hours.

.. _figTempHour9:

.. figure:: ./figTempHour9.png
   :align: center

   The simulation after nine hours.

.. _figTempHour13:

.. figure:: ./figTempHour13.png
   :align: center

   The simulation after thirteen hours.

.. _figTempHour16:

.. figure:: ./figTempHour16.png
   :align: center

   The simulation after sixteen hours.

.. _figTempHour21:

.. figure:: ./figTempHour21.png
   :align: center

   The simulation after twenty one hours.

Proposals for Topics of Projects
--------------------------------

1. geothermal heat pumps

   A geothermal heat pump compensates for the seasonal
   temperature fluctuations, to heat a home during the winter
   with the warmth of the earth stored during summer.

  * Use the Helmholtz phasor method to explain how deep
    the system must be buried.

  * Do a cost benefit analysis of installing
    a geothermal heat pump, comparing to heating
    a home with a furnace.

2. thickness of walls for isolation

   Native Americans of the southwest built homes with adobe walls
   of a critical thickness.  The walls delayed the flux of the sun
   until it was needed to heat the interior during the cold nights.

   * Model the thickness of the adobe walls.

   * Use actual material properties to select the optimal
     thickness for a 24-hour repetitive cycle.

   * Compare the results of your model against actual constructions.

3. study hibernating worms

   Birds migrate to survive the winter.  
   Worms have no wings but they do survive each winter.

   How deep do hibernating worms go to survive the winter?

   * Gather biological data on worms.

   * How deep does one have to dig to find worms?

   * Use the Helmholtz method to answer the question.

Exercises
---------

1. Consider the first proposed solution
   :math:`U(x, t) = X(x) \sin(\omega t - \theta(x))`.
   Use ``SymPy`` to substitute :math:`U(x, t)` in the heat equation.
   Argue why that method does not work.

2. Apply the Helmholtz phasor method to solve

   .. math::

      \frac{\partial u}{\partial t}
      = \frac{\partial^2 u}{\partial x^2}, \quad x \in [0, 1],

   subject to the boundary conditions:

   * :math:`u(0,t) = \sin(t)`,

   * :math:`u(1,t) = 0`.

3. We derived

   .. math::

      T(x, t) = T_0 + A e^{-x/\delta}
      \cos\left(\omega t - \frac{x}{\delta} \right).

   Show that by the substitutions
   :math:`(T(x,t) - T_0)/A = \widehat{T}` and
   :math:`\omega t = \widehat{t}`, 
   the model becomes

   .. math::

      \widehat{T}(x, \widehat{t}) = e^{-x/\delta}
      \cos\left(\widehat{t} - \frac{x}{\delta} \right),

   and thus dependent only on the penetration depth :math:`\delta`.

Bibliography
------------

1. Charles R. MacCluer, section 11.5 of
   *Industrial Mathematics.  Modeling in Industry, Science, and Government*,
   Prentice Hall, 2000.