Basically my question is why, when we go from classical solutions to weak solutions, we have to use energy inequalities instead of equalities. More preciscely, consider the Navier-Stokes equations
\begin{equation} \partial _t \rho + \text{div}(\rho u) = 0 \ \ \ \text{in } I \times \Omega \end{equation}
\begin{equation} \partial _t (\rho u) + \text{div}(\rho u \otimes u) + \nabla p(\rho) - \text{div}S(\nabla u) = \rho f \ \ \ \text{in } I\times\Omega \end{equation}
with
\begin{equation} u = 0 \ \ \ \text{on } \partial \Omega, \end{equation}
\begin{equation} \text{div}S(\nabla u) = \mu \Delta u + (\lambda + \mu)\nabla \text{div}u \end{equation}
For classical solutions, the momentum equation is multiplied by $u$ and by the continuity equation we obtain (using integration by parts)
\begin{equation} \int _{\Omega} \left(\frac{1}{2}\rho |u|^2 + P(\rho)\right)(t)dx + \int_I\int_{\Omega} \mu |\nabla u|^2 + (\lambda + \mu)|\text{div}u|^2dxdt = \int _{\Omega} \left(\frac{1}{2}\rho |u|^2 + P(\rho)\right)(0)dx + \int _I \int _\Omega \rho f \cdot u dxdt \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \, (*) \end{equation}
with $P(\rho) = \rho \int_1^\rho \frac{p(z)}{z^2}dz.$ Now for solutions of the weak formulation of the momentum equation, i.e.
\begin{equation} \int_I \int _\Omega \rho u \partial _t \phi + \rho u \otimes u : \nabla \phi + p(\rho)\text{div}\phi dxdt = \int_I \int_\Omega \mu \nabla u :\nabla \phi + (\lambda + \mu) \text{div}u \text{div}\phi - \rho f \phi dxdt \end{equation}
for all $\phi \in D(I\times \Omega)$ it says we have $(*)$ with $=$ replaced by $\leq$ but I don't see why. It says, that the equality is lost because of the weak lower semicontinuity in the term
\begin{equation} \int_I\int_{\Omega} \mu |\nabla u|^2 + (\lambda + \mu)|\text{div}u|^2dxdt, \end{equation}
so probably we have to use an approximation $(u_n)_n\subset D(I\times \Omega)$ to the weak solution $u$ as a test function and then pass to the limit, but this would only give us
\begin{equation} \int_I \int_\Omega \mu \nabla u :\nabla u_n + (\lambda + \mu) \text{div}u \text{div}u_n \end{equation}
and to apply weak lower semicontinuity we would rather need something like
\begin{equation} \lim_{n \rightarrow \infty}\int_I \int_\Omega \mu \nabla u_n :\nabla u_n + (\lambda + \mu) \text{div}u_n \text{div}u_n \geq \int_I\int_{\Omega} \mu |\nabla u|^2 + (\lambda + \mu)|\text{div}u|^2dxdt. \end{equation}
Thank you in advance for any hint on how to see this.