Spray combustion finds a wide range of application in gas turbines, internal combustion engines, industrial furnaces, etc. In turbulent spray combustion, liquid fuel is injected into the combustion chamber in the form of droplets. In order to improve the combustion efficiency and to reduce the thermal NOx emissions released during combustion, spray combustion could be operated in flameless mode. In a flameless combustion, the oxidizer is mixed with recirculated hot combustion gases to preheat and to dilute it. The dilution of oxidizer results in lower peak combustion temperature which reduces the NOx emissions and oxidizer preheating improves the thermal efficiency of combustion. In order to fully understand the combustion mechanics for its effective implementation in various applications, numerical simulations of flameless turbulent spray combustion are potentially useful because simulations are cost effective and serve as basis for further experimental studies based on the validation of numerical models. Turbulent spray combustion is a complex phenomenon involving two phases namely the gaseous phase and liquid phase. These two phases interact with each other through mass, momentum and energy transfer between them. This is complicated further by the interaction between the turbulence in the flow field and chemistry of the reacting species. Hence simplified models are necessary to simulate and understand the phenomenon of turbulent spray combustion. In this thesis, numerical validation study of turbulent ethanol spray flame using open-source software package OpenFOAM is carried out for the experiments done in Delft Spray-in-Hot-Coflow (DSHC) burner operated in flameless mode. The modeling approach used is Reynolds Averaged Navier Stokes simulations (RANS) with Eulerian-Lagrangian framework for the continuous phase and discrete phase respectively. Models like evaporation and turbulence models used in the sprayFoam solver are optimized for the spray combustion and validated with experimental data for one flame. The evaporation models studied are Gradient diffusion model and Spalding model. It is found that Gradient diffusion model gives better prediction of droplet properties at higher axial locations than Spalding model. The standard and realizable k-? model turbulence models comparative analysis showed that standard k-? model has much better gas phase temperature prediction than realizable k-? model due to the dependence of combustion model (Partially Stirred Reactor model) on the turbulence mixing frequency, ?/k. These optimized models are extended to simulation of HI and HIII flames. The peak gas phase temperature was under-predicted by the PaSR model. The results showed the importance of analyzing the different initial spray conditions.