Short question:
If the energy equation of an ideal plasma is written as follows:
$\begin{equation} \frac{\mathrm{d}p}{\mathrm{d}t}=-\Gamma p\nabla\cdot v-\left(\Gamma-1\right)\left( \nabla\cdot q-\Pi\left(\nabla v\right) \right), \end{equation}$
with p the pressure, $\Gamma$ the adiabatic index, v the velocity, q the heat flux through the boundary, and $\Pi$ the viscous stress tensor ($\Pi\left(\nabla v\right)$ the volumetric viscous heating rate), how would this look like in an expanding universe (preferably in a 3+1 notation)?
$\\$
$\\$
Long derivation and explanation of above formula:
If we are in the Lagrangian frame (a volume element co-moving with the fluid), the energy flow across the surface $-\nabla\cdot q$ plus the viscous heating rate $\Pi\left(\nabla v\right)$ must equal the heat input rate per unit volume $\rho\frac{\mathrm{d}Q}{\mathrm{d}t}$: \begin{equation} \rho\frac{\mathrm{d}Q}{\mathrm{d}t}=-\nabla\cdot q+\Pi\left(\nabla v\right), \end{equation} with $Q$ the heat unit per mass and $q$ the heatflux through the boundary.
The viscous heating rate $\Pi\left(\nabla v\right)$ comes from the viscous stress tensor $\Pi$ which is of the form \begin{equation} \Pi= \begin{pmatrix} 0 & a & b \\ a & 0 & c\\ b & c & 0 \end{pmatrix}, \end{equation} with, $a$, $b$, and $c$ scalars.
Lets assume there is no heat loss through the surface, but volume elements can interchange heat, so that the first law of thermodynamics \begin{equation} \mathrm{d}Q=p\mathrm{d}\left(\frac{1}{\rho}\right)+\mathrm{d}e, \end{equation} with $p\mathrm{d}\left(\frac{1}{\rho}\right)$ the PV-work per unit mass and $de$ the change in energy per unit mass can be substituted in above equation. We find \begin{equation} p\rho\frac{\mathrm{d}\left(\frac{1}{\rho}\right)}{\mathrm{d}t} +\rho\frac{\mathrm{d}e}{\mathrm{d}t}=-\nabla\cdot q+\Pi\left(\nabla v\right). \end{equation}
The factor $\frac{\mathrm{d}\left(\frac{1}{\rho}\right)}{\mathrm{d}t}$ can be rewritten: \begin{align} \frac{\mathrm{d}\left(\frac{1}{\rho}\right)}{\mathrm{d}t}&=\frac{-1}{\rho^2}\frac{\mathrm{d}\rho}{\mathrm{d}t}\\ &=\frac{-1}{\rho^2}\left(-\rho\nabla\cdot v\right)\\ &=\frac{1}{\rho}\nabla\cdot v. \end{align} With this we can express the rate of change of energy per unit volume as \begin{equation} \rho\frac{\mathrm{d}e}{\mathrm{d}t}=-p\nabla\cdot v-\nabla\cdot q+\Pi\left(\nabla v\right). \end{equation} Again, we can use the continuity equation, but this time to rewrite \begin{equation} \rho\frac{\mathrm{d}e}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\left(\rho e\right)+\rho e\nabla\cdot v, \end{equation} which leads to \begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\left(\rho e\right)=-\rho e\nabla\cdot v-p\nabla\cdot v-\nabla\cdot q+\Pi\left(\nabla v\right). \end{equation}
For an ideal gas we have \begin{equation} \rho e= \frac{p}{\Gamma-1}, \end{equation}
which leads to the equation on the top.