The Eulerian–Eulerian two-fluid model is widely used for computational fluid dynamics simulations of gas–solid flows. For non-isothermal flows, the averaged conservation equations solved in the two-fluid model require closures for drag, gas–solid heat transfer, pseudo-turbulent velocity fluctuations and the pseudo-turbulent heat flux (PTHF). However, the pseudo-turbulence terms are usually neglected in two-fluid simulations due to the lack of accurate correlations. With the increase in computational power, closures for these terms are now available from particle-resolved direct-numerical simulation (PR-DNS). Here, the PTHF closure as well as the heat-transfer closure (i.e., the Nusselt number) extracted from PR-DNS are implemented in the two-fluid thermal energy equation in OpenFOAM. The implementation is validated by comparing the simulation results with the PR-DNS data for the temperature profiles. Based on the analysis of the thermal energy budget, the PTHF can have a significant contribution and neglecting it can lead to large errors.