Interacting multi-fluids in FLRW cosmology? In the following, I only consider fully isotropic and homogeneous cosmology.
Because of local isotropy and homogeneity, the Robertson-Walker metric (and its two FLRW equations) imply a diagonal energy-momentum tensor, of a "perfect fluid" type.  The fluid is then described by two local quantities only: total energy density $\rho$ and total isotropic pressure $p$.  For a mixture of fluids, we could decompose these as sums:
\begin{align}
\rho &= \sum_i \rho_i,
& p &= \sum_i p_i. \tag{1}
\end{align}
The fluids must obey local conservation of energy (this comes from the Bianchi identity, so this equation can't be modified if we stay with full isotropy and homogeneity):
$$\tag{2}
\dot{\rho} + 3 \frac{\dot{a}}{a} (\rho + p) = 0.
$$
Now, for each fluid, we usually assume a linear macroscopic equation of state, where $w_i$ is generally assumed to be a constant:
$$\tag{3}
p_i = w_i \rho_i.
$$
The previous is standard FLRW cosmology.  There is no interaction between the various fluids into play.  Now I would like to go a bit beyond that to include fluid interactions, while keeping the full isotropy/homogeneity.  Lets talk about a mixture of two fluids to simplify things.  According to some papers on arXiv (for example : https://arxiv.org/abs/0803.1086), since $\rho = \rho_1 + \rho_2$ and $p = p_1 + p_2$, equation (2) could be splited into two parts:
\begin{align}
\dot{\rho_1} + 3 \frac{\dot{a}}{a} (\rho_1 + p_1) &= +\, Q, \tag{4a} \\[2ex]
\dot{\rho_2} + 3 \frac{\dot{a}}{a} (\rho_2 + p_2) &= -\, Q, \tag{4b}
\end{align}
where $Q$ is an arbitrary function.  Non-interacting fluids are associated to the special case $Q = 0$. Note that multiplying any of the equs (4) by $a^3 dt$ gives the following:
$$\tag{5}
dE = Q \, a^3 dt - p \, dV,
$$
where $E = \rho_1 \, a^3 = \rho_1 V$ is the fluid's energy in the volume $V = a^3$ (or $V \propto a^3$, if you prefer to include some constant, depending on the geometry of space).
So my question is this:
Can we interpret $Q \, a^3 dt$ in the right part of (5) as the heat $dQ$ exchanged between both fluids, in the time intervall $dt$?  If so, can we state that $dQ = T_1 dS_1 = T_2 dS_2$ for both fluids separately?  Or is that term better described as another interaction work $dW$ term?
 A: We should be a bit precise with terminology here. A perfect fluid is a fluid with stress tensor $T^{\mu \nu} = p g^{\mu \nu} + (\rho + p) u^\mu u^\nu$ (using -+++ signature) where $u^\mu$ is the fluid velocity vector satisfying $u_\mu u^\mu = -1$. Further, $\rho$ and $p$ are not independent; they are related by an equation of state. In full generality this would be of the form $F(\rho,p) = 0$ for some function $F$ (with some conditions imposed to exclude degenerate cases), but we can use the implicit function theorem to convert this to an equation of the form $p = G(\rho)$ where $G$ is now a function of a single variable. As a special case of this, we can restrict to linear equations of state, which are then $p = w \rho$ where $w$ is a constant number (called the equation of state parameter, though in practice people often abuse terminology and call $w$ the equation of state itself).
In cosmology, people often abuse terminology and call a fluid which has $\rho$ and $p$ given as functions of $a$ a perfect fluid. I find this confusing because such a fluid does not necessarily have an equation of state $p = g(\rho)$ except in the context of a particular background. For instance, if the fluid is sensitive to compression, it may have a bulk viscosity, which would add an effective pressure that depends on the background. Such a fluid is not actually perfect (though in some contexts this is not important).
Once you know the equation of state, you can almost define a temperature and entropy density for the fluid from the fundamental thermodynamic relation. Of course you can not do this exactly; for example, a conformal fluid has $p=\rho/3$ but to relate $p$ or $\rho$ to $s$ and $T$ you need to know the number of degrees of freedom $g$ of the fluid, with $\rho \propto g T^4$. This basically amounts to choosing an integration constant for the fundamental thermodynamic relation. Note that if your fluid has conserved quantities like particle number you'll need to account for those as well in the thermodynamic relation you write.
So far there is no cosmology here, only relativistic fluid dynamics. To do this in cosmology (at the level of a homogeneous, isotropic fluid) means to solve the Friedmann Equations, with whatever combination of fluids you like. In the absence of interactions between the fluids, this can be done and it isn't too hard either, but it's equivalent to just solving the equations with thermodynamics. More precisely, the comoving entropy $sa^3$ of a perfect fluid will be conserved; that is to say, the fluid cools adiabatically. So it is easy to find the entropy density as a function of the scale factor, and you can use that to find the energy density and pressure.
If it is surprising that this works, it probably shouldn't be: a perfect fluid is a fluid which equilibrates immediately, so cosmological evolution cannot push it out of thermal equilibrium. The microphysics rules the dynamics here and gravity merely provides it an expanding box to cool. The scale factor $a$ as a function of time will be whatever it needs to be to satisfy the Friedmann equations, and the fluid does not at all care because it will never be driven out of equilibrium. (If your fluid doesn't satisfy this, it isn't a perfect fluid.) There is nothing stopping you from putting as many non-interacting fluids as you like on the same FRW background; they will evolve independently of each other, each conserving its own entropy $s_i a^3$. The only point at which you need to add them together is to write down the Friedmann equations to solve for the scale factor as a function of time. If you really hate thermodynamics you can go through the exercise of writing down the Friedmann equations for your fluid with its complicated equation of state and solving them (even some textbooks do this), but you won't find anything new and will be forced to integrate everything numerically for any sufficiently nontrivial equation of state.
If you want to start allowing interactions between perfect fluids, then you have to actually do serious work (as if you considered an imperfect fluid which is homogeneous and isotropic). The reason is that the interactions you want to write down generally won't equilibrate quickly enough (relative to the Hubble scale) to be negligible. If they did, you would have fluids so tightly coupled that you would just treat them as a single fluid on the macroscopic scale. Specifically, to treat this, you want to write down the Boltzmann equations for the fluids. The equation you wrote down is a very simple version of this; in full generality with $n$ fluids the Boltzmann equations will look like
$$ \frac{\partial f_i}{\partial t} + L_i(f_i) = C_i (f_1, \ldots, f_n)$$
Here $f_i$ is the phase space distribution function for fluid $i$, $L_i$ is the Liouville operator, and $C_i$ are collision functionals that parametrize the interactions between the fluid. Of course, in the noninteracting case you have $C_i = 0$ and so each fluid evolves independently to follow Liouville's theorem.
In general these equations are very complex and only solvable numerically; the degree to which you are allowed to simplify it in your case depends on the nature of the collision functionals you want to write down. At the level of homogeneous, isotropic fluids, the phase space density is just a number and the only possibly collision functionals exchange energy (really entropy) between the various species of fluid, so they will look like what you wrote down (with $n$ fluids you will have $n$ functions on the right hand sides of each equation which sum to $0$). If you want to allow anisotropies in the fluid then you will need to consider more complicated collision functionals.
While we can write all this down, program it into a computer and get some numerical results, there are certainly some conditions that need to be satisfied for physical reasons. For sure, there are conditions on what sorts of equations of state $G$ are "physically admissible". The Einstein equations (which here simplify to the Friedmann equations) themselves do not care about the underlying theory that produces the stress tensor $T^{\mu \nu}$ so they can't tell us what is and isn't a valid perfect fluid. However, there are for instance the various energy conditions that arise as attempts to constrain what is a physically reasonable fluid. There is also the possibility of the Einstein equations developing an instability. There is a lot of literature here and still some open questions, but it's surely something that would require a more precise question to address in any real capacity.
