This work numerically investigates the turbulent partially premixed flames, which are ubiquitous in combustion powered devices. This combustion mode involves many physical complexities such as flame propagation in unevenly premixed mixture of fuel and oxidiser, turbulence/flame interaction in presence of mixture fraction gradients, triple flame configuration, etc. The fundamental mechanisms of these physical processes are yet to be further understood, posing significant modelling challenges. These issues are addressed in this thesis using a presumed joint probability density function (PDF) approach with both Reynolds Averaged Navier Stokes (RANS) and Large Eddy Simulation (LES) methodologies. This joint PDF is described by a parameter, mixture fraction, describing mixing and a reaction parameter, progress variable. The laminar flamelet concept is adopted to decouple chemistry and turbulence calculations for high computational efficiency. This modelling framework is validated using two experimental test cases in this study including a canonical lifted jet flame and a practical swirling flame, both exhibiting strong partial premixing features. The simulation results obtained for these validation cases show a robust model performance for a broad range of flow and mixing conditions with an attractive computational cost for practical interests. For the lifted jet flame case, two-dimensional (2D) axisymmetric steady RANS approach is used to compute the flame lift-off height showing very good agreement with the experimental measurements for a range of jet velocities and air-dilution levels. However, a substantial difference is found between the 2D unsteady RANS (URANS) results and experimental data for the flame transient evolution from its initial ignition to final stabilisation. The comparison for this transient evolution is improved significantly in the 3D URANS simulations suggesting that the third physical dimension, the azimuthal direction, plays an important role during the flame transient evolution. The following LES study further shows that the flame most-leading point appears to be in different azimuthal positions exhibiting a spiral -like trajectory as the flame propagates upstream towards the final lift-off height. The temporal variation of this leading point during this process is captured very well by the LES model in comparison with the experimental data. The validity of this LES model is then further assessed for a confined swirling flame with practical flow conditions. The simulation results are compared against an extensive experimental dataset including velocity, mixture fraction, temperature and major species mass fraction measurements, showing an overall good agreement at various locations inside the combustion chamber. The intermediate species mass fractions are also predicted reasonably well by the LES model despite that a general over-estimation is observed. Moreover, it is found that in the transport equation for the SGS variance of progress variable, the reaction and dissipation terms are predominant, substantially greater than the terms representing turbulent diffusion and production processes. Furthermore, the statistical correlation between mixture fraction and progress variable is found to be important for the RANS closure. However, this correlation is observed to be negligible for the sub-grid scale (SGS) of LES showing reasonably good model predictions for both simulated experimental configurations.