Nearly a year ago, I gave in the comments a non-answer that would justify calling the quantity
$$
\int u(x,t)^2~dx = E(t)
$$
the energy of the rod. Namely, one can think of minimizing the energy functional as satisfying a conservation law, and thus minimizing the energy leads you to solutions of the associated partial differential equations.
One year out, I think I have a more concrete, dimensional-analysis based answer. Let us examine the equation
$$
\partial_t u = -\alpha\partial_{xx}u
$$
for its dimensions. $u = u(x,t)$ represents temperature at the spacetime coordinate $(x,t)$, so it has units of temperature $T$. $\alpha$ here is the thermal diffusivity, which has units of length squared over time, $L^2/\tau$. Finally, not in the equation but lurking in the background is the kinetic energy, which is related to temperature by the equation
$$
E_k = \frac{3}{2}kT,
$$
where $k$ is the Boltzmann constant and has units of energy over temperature, $E/T$.
Let me use $[unit]$ to represent the dimensions of any units in any expression. So for example, the units in the heat equation check out: $\partial_t u$ has units of temperature over time, while $\partial_xx u$ has units of temperature over length squared, which gives us
$$
[T\tau^{-1}] = [L^2\tau^{-1}][TL^{-2}].
$$
Now, in our scenario the total length of the rod is fixed, so without loss of generality we may treat it as a constant. Dimensionally, this is saying we will take $[L] = [1]$, so it drops out of any of our equations. This is all fine, as we see in the heat equation the units of length cancel each other anyway.
So let us turn to the energy functional. In terms of units,
$$
\int u(x,t)^2~dx = \int [T^2]~d[L] = [T^2L] = [T^2].
$$
So the energy functional has units $[T^2]$, which isn't where we want to be (we want units of energy). However, the primary thing we do with taking the energy functional is minimize it. Minimizing the energy functional $E(t)$ is equivalent to minimizing the square root of the energy functional, $\sqrt{E(t)}$, which has units of $[T]$.
Still not there. But aha! Energy and temperature, $[E]$ and $[T]$, are related by a proportionality law! That is,the equation
$$
E_k = \frac{3}{2}kT
$$
tells us precisely how to convert from temperature to energy, and the conversion preserves order. So if we can minimize the square root of the "energy functional" $\sqrt{E(t)}$, which has units of $T$, then we automatically know how to minimize the actual energy, which according to our dimensional analysis must be something like
$$
\sqrt{\int \frac{9}{4}k^2 u(x,t)^2~dx } = \frac{3k}{2}\sqrt{E(t)},
$$
where the right-hand side is actually integrating the physical notion of energy squared. And now, the units do indeed check out: these are quantities with units of energy. There are surely other ways to define a natural notion of a physical energy functional using an integral, but the units of this definition work and it has nice mathematical properties (pretty much exactly those of $\sqrt{E(t)}$, which we know has an excellent mathematical theory). This revolves around the observation that since $k$ is a constant, minimizing the "energy" functional defined using units of heat is equivalent to minimizing what you would expect to be the "physical energy" functional.
Morally, the fact that $k$ is a constant (despite having units of $[ET^{-1}]$) gives us a temperature-energy equivalence law, which tells us that for dimensional analysis purposes $[E]$ and $[T]$ are indistinguishable. This is what the physicist in me needs to say to satisfy himself. The mathematician in me chooses units for temperature and energy so that $k = 2/3$, and then I just have the simple equivalence $E_k = T$, and now I happily exchange temperature and energy at will. And presumably this is what mathematicians of the past have done.
Treating $L$ like a constant also needs some justification, but I think if we think of $L$ as very small and argue that computing the energy of large homogeneous bodies consists of computing the energy on small pieces and summing, then we can justify that assumption as well.
In general, for any $n$ you can look for solutions involving terms where $(\text{degree in }x) + (\text{degree in }y) + 2 (\text{degree in }t) = n$.
To find such a solution, start with a polynomial $p_0(x,y)$ homogeneous of degree $n$. The solution will be
$$ \sum_{m=0}^{\lfloor n/2 \rfloor} \frac{t^m \Delta^m p_0}{m!}$$
For example, starting with $p_0(x,y) = x^2 y^3$ homogeneous of degree $5$
we get the solution ${x}^{2}{y}^{3}+t \left( 6\,{x}^{2}y+2\,{y}^{3} \right) +12\,{t}^{2}y$.
Best Answer
There are some physical limits which are often hidden in the usual presentations. The "physically reasonable" boundary condition is the heat flux boundary condition: there is a thermal conductivity coefficient $k$ for the boundary, and the system is in thermal contact with a temperature $u_{ext}$ at each point of the boundary. Then the boundary condition says that the outward heat flux at a point $x$ on the boundary is $k(u(t,x)-u_{ext}(t,x))$. This outward heat flux is $\nabla_x u \cdot n$. So you wind up with
$$\begin{cases} u_t(t,x) = \Delta u(t,x), & x \in \Omega,t>0 \\ (\nabla_x u \cdot n)(t,x) = k(t,x) (u(t,x)-u_{ext}(t,x)), & x \in \partial \Omega \\ u(0,x)=f(x). & \end{cases}.$$
The Dirichlet boundary condition occurs when you send $k \to \infty$; in this case the boundary is in perfect thermal contact with the heat reservoir, and so $u$ is forced to be equal to $u_{ext}$. The homogeneous Neumann boundary condition occurs when you send $k \to 0$; in this case the boundary is ideally isolated from the heat reservoir, and so the total heat in the system is held fixed. Note also that when we take material properties into account, $k$ could be a source of nonlinearity, in that the heat flux could be a more complicated function of $u$ and $u_{ext}$.
Responding to some of your middle remarks: