Considering the future coming more stringent emission limits, the working condition of the combustors will be more rigorous, and the problems will be more serious. To meet this high industry demand, a deeper and better understanding of the physics of the turbulent combustion process is wanted. In doing this, more accurate computational studies are necessary. Eventhough, direct numerical simulations (DNS) of Navier-Stokes equations can predict turbulence without models; these types of simulations are restricted to very simple geometries. On the otherhand the large eddy simulations (LES) are now viewed as a promising tool to address turbulent combustion flows where classical Reynolds-averaging numerical simulations (RANS) approaches have proved to lack precision or where the intrinsically unsteady nature of the flow makes RANS clearly inadequate. In this thesis large eddy simulation techniques are considered. This work consists of two parts. First, this thesis studies the subgrid-scale (SGS) modelling for large eddy simulations. The aims are to gain deeper understanding of the sensitiveness of the SGS modelling and predict the simple and better model for the further studies. For this purpose, effect of different filter width formulas on Smagorinsky model, robustness of the dynamic model, behaviour of different dynamic procedures, behaviour of mixed models and one equation models are studied. Nevertheless all these models are not capable of predicting the anisotropic nature of the flow phenomena. In this work anisotropic one equation model based on the subgrid scale (SGS) kinetic energy is developed. This model considers an anisotropic eddy viscosity formulation. All these models are tested by considering a turbulent channel flow. Finally some models are tested on rotating channel flow. It is shown that the proposed model achieves a good agreement with the experimental data and significantly outperforms isotropic models. Nevertheless Lagrangian dynamic model also predicts good flow phenomena. The Second part consists of the large eddy simulation of turbulent premixed combustion processes by considering the flame surface density and artificial thickened flame model. For this purpose a simple technique for the density coupling is proposed for the extension of the incompressible solver to variable solver. In this technique, first divergence of the velocity is calculated and then density is updated by using this value. Divergence of the velocity is also used for the treatment of convective outlet boundary condition. Handling of the flame wall interaction is presented. Effect of different flame wrinkling models is studied. For simulations bluff-body stabilized flame is considered. Comparisons between the experimental and numerical results have showed that the numerical method implemented is able to provide promising agreement.