In fluid dynamics, especially in modeling, there are different flavors for including heat.
First of all, heat can be a passive tracer which does not influence the flow, this basically means, that you solve the scalar transport-diffusion equation for the temperature (or energy if you which), which reads, in your notation
∂T∂t+(u⋅∇)T=α∇2T+Q
where α is thermal diffusion coefficient and Q are heat sources (which may be dependent on T for radiation).
The density of fluids is temperature-dependent. You should incorporate this in the Navier-Stokes equations. If the density differences are only relevant for the buoyant terms (e.g. gravity), it is enough to include the effect in f, by the Boussinesq approximation. This basically means that you linearize the ρ(T) curve and include the force term as βg(T−Tref) with β a thermal expansion coefficient.
The next step in complexity, arises when ρ changes such, that it takes over control for more terms. The Navier-Stokes equation you gave, assumes constant density, which therefor drops out of the equations. This is generally not true, and you have to take into account the variable density in all terms.
This post imported from StackExchange Physics at 2014-03-17 04:25 (UCT), posted by SE-user Bernhard