The thesis is focused on the current state-of-the-art of modeling turbulent single-phase and turbulent two-phase stratified flows. This report is meant to identify and validate the most suitable candidate model for inclusion in the multiphase flow models that are being developed in the Scientific Computing group of the Delft Institute for Applied Mathematics. For that reason, as an initial step the results of single phase turbulent pipe flows which are simulated using DNS and LES methods are compared with the results of Eggels. In two-phase flows, if the gas phase is turbulent the liquid phase will become turbulent as well. If the transition occurs from stratified flow to stratified-wavy flow, the interfacial momentum transfer varies due to the existence of waves at the interface. This process makes modeling of the momentum transfer complicated. In general, when the effects of surface tension are negligible the equations for the two-phase flow and for the single-phase flow are identical, the only differences between two-phase and single phase flow equations are the variable density and viscosity. Therefore, the influence of the interface and the momentum transfer between both phases can be ignored and a simple single-phase flow model combined with an interface model can be used as an initial approximation while concentrating on modeling the turbulence in both phases away from the interface. For this reason, the one-fluid model needs to be introduced (refer to Appendix \ref{App:s62}) in order to obtain results for multiphase flows using classical single-phase flow models. In this project, two different types of numerical techniques, namely DNS and LES are chosen to estimate the computational resources of a turbulent single-phase pipe flow test case with a friction Reynolds number of $Re_*=360$. Accordingly, computational complexities of different techniques are analyzed in detail. The estimation procedure of the problem complexity (i.e., the required number of total grid cells) for turbulent single-phase flows gives an underestimate for the number of unknowns of turbulent two-phase flows. The comparison of computational costs showed that Direct Numerical Simulation (DNS) is possible for turbulent two-phase stratified pipe flows only for low Reynolds numbers. For high Reynolds number flows, DNS is not feasible because the current implementation of the algorithm is not parallelized and the computational resources of the Scientific Computing group are limited. Because of these reasons, large eddy simulation (LES) is considered to be the promising technique as the computational resources required for DNS become excessive for higher Reynolds number the serial code. Therefore, LES needs to be investigated elaborately for turbulent two-phase stratified pipe flows in future. At this point, several numerical simulations are performed to take first steps towards simulation of turbulent two-phase flows using the one-fluid model. However, it is due to limitations in time, turbulent two-phase flow simulations are not performed and only turbulent single-phase flows are considered in this thesis. The numerical results for the Poiseuille flow are obtained both in Cartesian and cylindrical coordinates to verify the variable viscosity formulation before analyzing turbulent flows. The algorithm developed in the Scientific Computing group is improved with necessary periodic boundary conditions for the discretization of the equations that describe turbulent single-phase flow in a circular pipe geometry. First, these boundary conditions are introduced into the algorithm and then results of numerical simulations in 2D (channel flow) and 3D (axisymmetric pipe flow) are validated by comparing them with theoretical values in Section 5. Then, the variable viscosity formulation is incorporated into the algorithm to take a step towards LES computations. In addition to this, the subgrid scale (SGS) parametrization and the Smagorinsky model are utilized for treatment of SGS turbulence which constitutes the basis of LES. The present numerical results for both DNS and LES illustrate that they are in agreement with the results of Eggels. Moreover, both methods are capable of simulating the problem within a reasonable amount of time and accuracy. It is also shown that choosing a relatively smaller pipe length (because of the restrictions imposed by the serial code) than the one chosen by Eggels has no significant effect on the resulting velocity profiles.