The electroluminescence from organic light-emitting diodes can be predicted with molecular-scale resolution using three-dimensional kinetic Monte Carlo (3D KMC) simulations [M. Mesta et al., Nat. Mater. 12, 652 (2013)]. However, around and below the built-in voltage KMC simulations are computationally inefficient. 3D master-equation (3D ME) simulation methods, which are fastest for low voltages, are so far mainly available for describing unipolar charge transport. In such simulations, the charge-carrier interactions are treated within a mean-field approach. It is not a priori evident whether such simulations, when applied to bipolar devices, can be extended to include the Coulomb attraction between the individual electrons and holes, so that charge-carrier recombination is sufficiently well described. In this work, we develop a systematic method for extending 3D ME simulations to bipolar devices. The method is applied to devices containing materials with Gaussian energetic disorder, and validated by a comparison with the results of 3D KMC simulations. The comparison shows that the 3D nonuniformity of the molecular-site-resolved carrier concentration and the one-dimensional layer-averaged profile of the recombination rate are fully retained, and that the 3D nonuniformity of the molecular-site-resolved recombination rate is fairly well retained.