We formulate new multi-phase convective heat transfer equations by combining the three-dimensional (3D) Navier-Stokes equations, the energy equation and the Cahn-Hilliard equation for the phase field variable phi(X, t). The density, viscosity, heat capacity and conductivity are functions of phi(X, t). The equations are solved in time with a splitting scheme that decouples the flow and temperature variables, yielding time-independent coefficient matrices after discretization, which can be computed during preprocessing. Here, a spectral element method is employed for spatial discretization but any other Eulerian grid discretization scheme is also suitable. We test the new method in several 3D benchmark problems for convergence in time/space including a conjugate heat transfer problem and also for a realistic transient cooling of a 3D hot object in a cavity with a moving air-water interface. These applications demonstrate the efficiency of the new method in simulating 3D multi-phase convective heat transfer on stationary grids, different modes of heat transfer (e.g. convection/conduction), as well as its robustness in handling different fluids with large contrasts in physical properties. (C) 2014 Elsevier Ltd. All rights reserved.