Turbulent spray flames can be found in many applications, such as Diesel engines, rocket engines and power plants. The many practical applications are a motivation to investigate the physical phenomena occurring in turbulent spray flames in detail in order to be able to understand, predict and optimise them. Turbulent spray flames are a two-phase flow system consisting of a turbulent reacting continuum phase and a dispersed liquid phase contained in the continuum phase. In this thesis the interaction between these two phases is being studied, and specifically the effects of vaporisation and the effects of the presence of droplets on the turbulence characteristics of the continuum phase (the two-way coupling effect). The main goal of the research described in this thesis is to develop and evaluate stochastic Lagrangian models for turbulent spray flames that are characterised by significant or dominant two-way coupling effects. In the first part of the thesis, isolated droplet vaporisation models are studied and several models are evaluated. The problem of droplet interaction is briefly explained using some simple simulations and general theoretical ideas and the concept of group combustion is introduced. In the second part of the thesis the hybrid Lagrangian-Lagrangian method is presented. The continuum phase is described by the transport equation for the joint mass density function (MDF) of velocity and composition, which is solved indirectly using a Monte Carlo particle method where Langevin type equations are solved. This method is augmented by a finite volume method where the mean Eulerian transport equations for momentum, Reynolds stresses and turbulence dissipation are solved, making it a hybrid approach for the continuum phase. The dispersed phase is described solely by the droplet MDF which is also solved using a Monte Carlo method. A novel modification of the turbulence model for the pressure rate of strain in the Eulerian transport equation for the Reynolds stresses is presented that takes into account the effect of the presence of particles and the effect of mass transfer. The models in the Langevin equation in the Lagrangian formulation are then modified in order to be consistent with the p ressure rate of strain model in the Eulerian formulation. The modification of the diffusion term in the Langevin equation implies a modification of Kolmogorov's hypothesis regarding the first order Lagrangian structure function. An investigation of the transport equation for the scalar dissipation rate shows that Lagrangian mixing models are affected by liquid particle vaporisation. Other theoretical developments are the extension of the Generalised Gradient diffusion model and the Daly Harlow model for the triple correlations and an exact expression for the effect of mass transfer on the turbulence dissipation. In the third part of this thesis the results of numerical simulations of two test cases are presented. In the first test case an Eulerian-Lagrangian approach was used and the influence of the modification of the model for the pressure rate of strain by the presence of particles is investigated for an axially symmetric non-burning turbulent dispersed spray and compared with experimental data. Small improvements in turbulence quantities over the conventional turbulence and two-way coupling model are noticed when the novel modification of the pressure rate of strain is being used. In the second test case, simulations of an axially symmetric turbulent spray flame are performed to illustrate the performance of the complete hybrid method proposed in this thesis. A good agreement with experimental data was found using the hybrid method. In the mass transfer controlled spray flame studied in this thesis, taking into account the reduction of the drag coefficient and the heat and mass transfer coefficient due to mass transfer turned out to be crucial to the performance of the method.