The objective of the research presented in this thesis is development and validation of predictive models or modeling approaches of liquid fuel combustion (spray combustion) in hot-diluted environments, known as flameless combustion or MILD combustion. The goal is to combine good physical insight, appropriate numerical methods and good software development in the context of the general framework of Computational Fluid Dynamics (CFD). For model validation to be possible, availability of relevant, accurate and complete experimental datasets is important. For this study the database of Delft Spray-in-Hot-Coflow(DSHC) flames has been employed. Using different combinations of model components and approaches, the modeling has been developed and tested at three levels of complexity and accuracy. Depending on the chosen combination of turbulence and combustion models, the modeling methods used in this thesis are the following: 1) Reynolds Averaged Navier-Stokes (RANS) approach for turbulence modeling with Steady Flamelet Model (SFM) approach for combustion; 2) Transported Probability Density Function (TPDF) approach for turbulence modeling with Flamelet Generated Manifolds (FGM) for combustion; 3) Largy Eddy Simulation (LES) for turbulence with FGM for combustion. The CFD codes that are used for these methods are respectively the commercial software ANSYS Fluent (version 15.0), the in-house code “PDFD” and the open source software OpenFOAM. The first stage of the research aimed at exploring the characteristics of the target flames, and gaining understanding of the relevance of each model component in order to clarify the directions of further model improvements. Existing standard models (RANS/SFM) within the commercial software were employed for this purpose. The main outcome of this stage of the study is that the SFM fails in predicting the lifted-off phenomena of the DSHC flames, and it is concluded that a more sophisticated Turbulence-Chemistry Interaction (TCI) model is required. In the second stage of the study, advancements on both turbulence and combustion models have been made, respectively with TPDF and FGM. TPDF was used because it offers a more detailed description of statistical properties and FGM because it offers the possibility to include finite rate chemical effects in the frame of TPDF at affordable cost. In this stage of the study the focus was on simulations of the dilute spray region. Ignoring the dense region close to the atomizer, the simulation starts from a certain distance downstream, preferably a region where experimental data are available or good insight in local conditions is available by other means. Because of the careful specification of the boundary conditions of spray as well as gas phase, the flow and combustion development downstream were correctly reproduced by this approach, including the lift-off height. Next, based on this platform, parametric studies on many important aspects have been carried out, for example, the influence of the droplet initial temperature, a comparison between different model variants on evaporation etc. Good agreement with experimental data was achieved and useful information on model performance was obtained. However, also limitations of the modeling approach were identified. The most important one is that the use of an adiabatic approach disregarding the energy consumed in the evaporation of liquid (adiabatic FGM table) over-predicts the gas phase temperature in the region of intense droplets evaporation. In the third stage of this study we used Large Eddy Simulation of turbulence in combination with non-adiabatic FGM for evaporation and combustion. The subgrid scale fluctuations of the FGM scalars were modeled using assumed-shape PDFs. The CFD platform for doing so is the open source CFD package — OpenFOAM. Since the desired FGM model is not available in the public released version of OpenFOAM, the first step in this third stage was development and validation of a new implementation of FGM. A well-documented gaseous lifted flame in hot-diluted coflow was used as test case. In this study, it was found that the auto-ignition process is highly sensitive to the model constants that are used to evaluate the variances of mixture fraction and progress variable. Dynamic procedures have been developed for the determination of these constants in the context of LES. Using the developed dynamic models significant improvement of the prediction of the flame lift-off height has been achieved. The second step of LES/FGM method development was the extension of the implemented FGM method to spray combustion. Two major model developments were completed in this step — the Condition Droplet Injection Model (CDIM) and the non-adiabatic FGM method. The CDIM was proposed to take into account the influence of flash-boiling atomization of the DSHC flames. The comparison of CIDM with a conventional droplet injection model showed superior performance of the CDIM in modeling of the DSHC flames. Furthermore, with the addition of enthalpy deficit as an extra dimension of the FGM table, the prediction of the gas phase temperature was significantly improved. The third step of modeling stage three was to apply the developed LES/FGM method to a wide range of test cases including the hot and cold coflow cases to test the range of applicability of models, and to some virtual (i.e. not yet studied experimentally) cases in order to gain deeper insight on the involved phenomena. The simulations showed that the developed LES/FGM method is able to correctly capture major features of spray flames under different coflow conditions, for example the flame width and flame lift-off height. Independent parameter study of the coflow temperature and oxygen concentration clarified the effects of these two important factors, and suggested better operation conditions to achieve a strict MILD spray combustion condition. Simulation and comparison of the hot coflow (Hii) and cold coflow (Aii) cases revealed very interesting mechanisms that determine different flame structures and their transition in spray combustion. The “single” and “double” flame structures of these two cases have been successfully reproduced. Detailed analysis showed that the actual flame topologies have been over-simplified by their conventional names (single/double flame) used in literature, which are mainly based on the experimental observation. It was found that the appearance of multi-flame or multiple reaction regions in spray combustion are resulting from the disparity of time scales of the processes involved. And it was further demonstrated that by matching important time scales, similar flame structure can be achieved under considerably different operating conditions. In summary, three different levels of modeling approaches for spray combustion have been developed and tested. These approaches can be employed in different applications based on the consideration of the required accuracy and computational affordability. The RANS/SFM is the cheapest one in terms of computation cost, but the interpretation of the results should take into account that a number of important aspects such as lift-off height are not well predicted. The LES/FGM method provides most accurate results compared to other two, but is also most computationally demanding. The method using the last approach offers the most promises for a better understanding of the MILD spray combustion and design of clean and efficient combustion technology.