In this paper we consider, for modelling and simulation, a non-isothermal turbulent flow laden with non-evaporating spherical particles which exchange heat with the surrounding fluid and do not collide with each other during the course of their journey under the influence of the stochastic fluid drag force. In the modelling part of this study, a closed kinetic or probability density function (p.d.f.) equation is derived which describes the distribution of position x, velocity v, and temperature θ of the particles in the flow domain at time t. The p.d.f. equation represents the transport of the ensemble-average (denoted by 〈 〉) phase-space density 〈W(x, v, θ, t)〉. The process of ensemble averaging generates unknown terms, namely the phase-space diffusion current j = βv〈u′W〉 and the phase-space heat current h = βθ〈t′W〉, which pose closure problems in the kinetic equation. Here, u′ and t′ are the fluctuating parts of the velocity and temperature, respectively, of the fluid in the vicinity of the particle, and βv and βθ are inverse of the time constants for the particle velocity and temperature, respectively. The closure problems are first solved for the case of homogeneous turbulence with uniform mean velocity and temperature for the fluid phase by using Kraichnan’s Lagrangian history direct interaction (LHDI) approximation method and then the method is generalized to the case of inhomogeneous flows. Another method, which is due to Van Kampen, is used to solve the closure problems, resulting in a closed kinetic equation identical to the equation obtained by the LHDI method. Then, the closed equation is shown to be compatible with the transformation constraint that is proposed by extending the concept of random Galilean transformation invariance to non-isothermal flows and is referred to as the ‘extended random Galilean transformation’ (ERGT). The macroscopic equations for the particle phase describing the time evolution of statistical properties related to particle velocity and temperature are derived by taking various moments of the closed kinetic equation. These equations are in the form of transport equations in the Eulerian framework, and are computed for the case of two-phase homogeneous shear turbulent flows with uniform temperature gradients. The predictions are compared with the direct numerical simulation (DNS) data which are generated as another part of this study. The predictions for the particle phase require statistical properties of the fluid phase which are taken from the DNS data. In DNS, the continuity, Navier–Stokes, and energy equations are solved for homogeneous turbulent flows with uniform mean velocity and temperature gradients. For the mean velocity gradient along the x2- (cross-stream) axis, three different cases in which the mean temperature gradient is along the x1-, x2-, and x3-axes, respectively, are simulated. The statistical properties related to the particle phase are obtained by computing the velocity and temperature of a large number of particles along their Lagrangian trajectories and then averaging over these trajectories. The comparisons between the model predictions and DNS results show very encouraging agreement.