Creating a reliable and accurate Red Blood Cell (RBC) aggregation model for small and midsize arteries and veins is still an active research subject with more in focus with a multi-scale approach including mesoscale effects. Better understanding the RBC aggregation requires a multi-phase and multi-scale approach for simulating blood with Newtonian and non-Newtonian parts. In our proposed work, viscosity, shear rates, phase distributions and volume fractions with a range of hematocrit levels of RBC are calculated using the depletion interaction theory for two-phase blood flow simulation and compared with the numerical and experimental data in literature. In addition, thermal effects are modeled using energy equations and changes in RBC aggregation are studied with respect to thermal variations. Two-phase fluid-fluid model is used including inter-phase momentum exchange. A new shape factor is proposed for the coupling effects on drag and lift forces. Finally, total interaction energy of RBCs, hematocrit levels of blood at varying temperatures and effects of temperature on viscosity and relative apparent viscosity are computed at varying shear rates and compared with the existing data in literature.