The accurate prediction of turbulent mixing induced by Rayleigh–Taylor (R–T), Richtmyer–Meshkov (R–M) and Kelvin–Helmholtz (K–H) instabilities is very important in understanding natural phenomena and improving engineering applications. In applications, the prediction of mixing with the Reynolds-averaged Navier–Stokes (RANS) equation remains the most widely used method. The RANS method involves two aspects, i.e. physical modelling and model coefficients. Generally, the latter is determined empirically; thus, there is a lack of universality. In this paper, inspired by the well-known Reynolds decomposition, we propose a methodology to determine the model coefficients with the following three steps: (i) preset a set of analytical RANS solutions by fully using the knowledge of mixing evolutions; (ii) simplify the differential RANS equations to algebraic equations by imposing the preset solutions to RANS equations; (iii) solve the algebraic equations approximately to give the values of the entire model coefficients. The specific application of this methodology in the widely used K–L mixing model shows that, using the same set of model coefficients determined from the current methodology, the K–L model successfully predicts the mixing evolutions in terms of different physical quantities (e.g. temporal scalings and spatial profiles), density ratios and problems (e.g. R–T, R–M, K–H and reshocked R–M mixings). It is possible to extend this methodology to other turbulence models characterised with self-similar evolutions, such as K-$\epsilon$ mixing models.