The graduation project was conducted at the CFD department of Nuclear Research And Consultancy Group (NRG) in Petten. The modeling and simulation of Taylor bubble flow using CFD can contribute significantly to the topic of nuclear reactor safety and in particular, in the emergency cooling of nuclear reactors during a loss of coolant accident, or in the U-tubes of a stream generator during a pipe rupture. To achieve an accurate representation of the gas-liquid interface for high values of Reynolds number, a general interfacial turbulence model should be developed which adapts to local conditions automatically. Direct Numerical Simulation (DNS) of relevant large interface two-phase turbulence has the potential to contribute to this, as it can produce more refined insight while being complementary to experimental data. The current graduation project illustrates a simulation approach towards DNS of turbulent co/countercurrent Taylor bubble flow. It comprises a continuation of the novel simulation strategy indicated by (Frederix et al., 2020) in which LES of co-current turbulent Taylor bubble flow using OpenFOAM were indicated and the authors concluded that LES mesh resolution is not sufficient to capture accurately the gas-liquid interface, velocity fluctuations, and bubble disintegration rate. To counter this, in the current work, a Basilisk code is developed which due to its adaptive local grid refinement and its accurate solution of advection equation, comprises a better choice than OpenFOAM for DNS in two-phase flows (via the settings of (Hysing et al., 2009) for laminar bubble flow and (Shemer et al., 2007) for turbulent co-current flow) since it captures sharper the interface and reduces the computational cost significantly. Except for OpenFOAM, Basilisk is also successfully validated against ANSYS (Araujo et al., 2012) for the laminar Taylor bubble flow. Last but not least, the work is further extended to the simulation of laminar and turbulent counter-current Taylor bubble flows in which the effect of the choice of the pipe diameter and the initial bubble length on the bubble’s decay rate is analyzed for the experimental setting of (Mikuz et al., 2019). Overall, despite the lack of fully DNS quality, this study extends the work of (Frederix et al., 2020) and provides further insight into the performance of Basilisk in two-phase flows at a reasonable computational cost. The results and conclusions from the current study may contribute to the development of low-order turbulence models and the validation of more general two-phase modeling strategies.