Similar to single-phase flows, multiphase systems with a constant heat flux at the wall are common in practice. Unlike single-phase flows, however, the numerical implementation of a constant boundary flux is non-trivial due to the coupling between phases - i.e., the partition of the total flux between the phases can vary in space and time. A numerical technique for modeling such a boundary condition is proposed here and verified via simulations of gas-solid flows. While discrete-particle simulations of monodisperse particles are considered here, the technique can be extended to include radiation, polydisperse systems, and/or a continuum representation of the solids phase. (C) 2016 Elsevier B.V. All rights reserved.