We develop a two-fluid model (TFM) for heat transfer in dense non-Brownian suspensions. Specifically, we propose closure relations for the inter-phase heat transfer coefficient and the thermal diffusivity of the particle phase based on calibration against experimental data. The model is then employed to simulate non-isothermal flow in an annular Couette cell. We find that, when the shear rate is controlled by the rotation of the inner cylinder, both the shear and thermal gradients are responsible for particle migration. Within the TFM framework, we identify the origin and functional form of a "thermo-rheological" migration force that rationalizes our observations. Furthermore, we apply our model to flow in eccentric Couette cells. Our simulations reveal that the system's heat transfer coefficient is affected by both the classic shear-induced migration of particles and the newly identified thermo-rheological migration effect. Finally, we employed the proposed computational TFM framework to analyze electronics cooling by forced convection for microchannel cooling. We used a suspensions of high thermal conductivity (Boron Nitride) particles in a 3M Fluorinert FC-43 cooling fluid. Three-dimensional simulations were run to quantify the temperature distributions under uniform heating (5 W) and under hot-spot heating (2 W/cm^2) conditions. A 100 K junction level temperature improvement (enhanced thermal spreading) was seen for hot-spot heating and 15 K was observed for uniform heating, demonstrating the enhanced cooling capabilities of dense particulate suspensions of high-conductivity particles, over a clear FC-43 fluid.