Potential energy in continuum physics: Why is this $0$ in the surface? At the start of the continuum mechanics, we were calculating the potential energy. Here is the part that's not clear to me:
Let $V$ be the potential energy and $U$ be its density function. Then:
$$V=\int U \mathrm{d}V$$
We have that
$$\nabla V=-\underline{F}^e$$
$$\nabla U = -\underline{f}^e$$
$$\delta U=-\underline{f}^e \delta\underline u$$
Where $\underline F^e$ is the elastic force, $\underline f^e$ it's density and $\underline u$ is the displacement vector.
$$\delta V = -\int \underline f ^e \delta \underline u \mathrm{d}V$$
$$= -\int \partial_l\,\sigma_{kl}\,\delta u_k\, \mathrm{d}V$$
$$= -\int \left[\partial_l \, (\sigma_{kl} \,\delta u_k)-\sigma_{kl} \, \partial_l \, \delta u_k\right] \mathrm{d}V$$
$$= -\oint \sigma_{kl} \,\delta u_k \mathrm{d}A_l+\int \sigma_{kl} \, \partial_l \, \delta u_k \mathrm{d}V$$
Because $\sigma_{kl}=0$ on the surface:
$$=\int \sigma_{kl} \, \partial_l \, \delta u_k \mathrm{d}V$$ 
Q1: What does exactly the $\delta U$ and $\delta \underline u$ mean?
Q2: Why is the $\sigma_{kl}=0$ on the surface?
 A: The quantities $\delta U$ and $\delta V$ are the infinitesimal variations of the potential density and potential respectively. There is an alternative derivation (from Landau & Lifshitz) that one can use that may help your understanding with the second question.
Consider the work required to deform a volume $\Omega$ (enclosed by surface S) by some infinitesimal displacement $\delta \textbf{u}$, 
$$W = \iiint_\Omega \delta W d\Omega = \iiint_\Omega \textbf{b} \cdot \delta \textbf{u} d\Omega + \iint_S \left(\sigma : \hat{n} \right) \cdot \delta \textbf{u} dS.$$
Note we are just simply writing this as a sum of body and surface work. Here $\textbf{b}$ is the body force density and $\sigma$ is the stress tensor. The notation $\sigma : \hat{n}$ means to contract the tensor onto the surface normal ($\sigma_{ij} \hat{n}_j$). Using the divergence theorem, one finds
\begin{align}
W &= \iiint_\Omega \textbf{b} \cdot \delta \textbf{u} d\Omega + \iiint_\Omega \nabla \cdot \left[\sigma : \delta \textbf{u} \right] \cdot \delta \textbf{u} d\Omega \\ \nonumber
&= \iiint_\Omega \textbf{b} \cdot \delta \textbf{u} d\Omega + \iiint_\Omega \left(\nabla : \sigma \right) \cdot \delta \textbf{u} d\Omega +\iiint_\Omega \sigma : \delta \left(\nabla \cdot \textbf{u} \right) d\Omega\\\ \nonumber
&= \iiint_\Omega \left(\textbf{b} + \nabla : \sigma \right) \cdot \delta \textbf{u} +\iiint_\Omega \sigma : \delta \left(\nabla \cdot \textbf{u} \right) d\Omega.\\ \nonumber
\end{align}
Note immediately that the first term $\textbf{b} + \nabla : \sigma$ is zero since that is identically the condition of static mechanical equilibrium. This is also called Cauchy's first law of momentum or the stress divergence equation. So, therefore
$$W = \iiint_\Omega \sigma : \delta \left(\nabla \cdot \textbf{u} \right)d\Omega.$$
which is exactly what you have above. Relating this back to $W = \int_\Omega \delta W d\Omega$ shows that $\delta W = \sigma_{ij} \delta \varepsilon_{ij}$
The constraint $\sigma_{kl}|_S = 0 $ is sometimes called a stress-free boundary condition. In continuum mechanics, this is a natural surface condition which indicates you do not have any shearing forces because they would violate mechanical equilibrium. Note that
$$\sigma_{kl} = C_{klij} \varepsilon_{ij} = \frac{1}{2} C_{klij} \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}  \right)  = 0 $$
where $C_{ijkl}$ is the elastic stiffness tensor.
