Turbulent mixing induced by interfacial instabilities, such as Rayleigh–Taylor (RT), Richtmyer–Meshkov (RM), and Kelvin–Helmholtz (KH) instabilities, widely exist in natural phenomena and engineering applications. On the one hand, the Reynolds-averaged Navier–Stokes (RANS) method, mainly involving physical model and model coefficients, is still the most viable approach in application. On the other hand, predicting different mixing problems with the same physical model and model coefficients—defined as “unified prediction” in this paper—is the basis for practice because (1) different instabilities usually exist simultaneously in a flow system and are coupled to each other; (2) mixing processes involve a wide range of parameters (e.g., time-dependent density ratio and acceleration history, etc.). However, few models can achieve such a unified prediction. Recently, we proposed a RANS route to realize this unified prediction by setting model coefficients to match the given physical model. This study attempts to apply this to the widely used BHR2 model to achieve unified predictions of different turbulent mixing problems, including basic problems (i.e., classical RT, RM, and KH mixing) and complex problems (i.e., re-shocked RM, tilted-RT, and spherical implosion mixing). Good agreement between experiments, large-eddy simulations, and RANS results were obtained. The temporal evolution of mixing width and spatial profiles of important physical quantities are presented. Based on our achievements of the k – L and k−ε models for unified predictions, the success of BHR2 model further confirms that our RANS route is robust for different turbulent mixing models and may be expanded to other fields.