The major problem regarding numerical simulation of multiphase flows in microfluidic systems are discrepancies in the results with the ones obtained by the same numerical approaches in larger scales. Several potential causes such as the effect of confining boundaries, large curvature of interfaces or small sizes of dispersed phases are reported elsewhere. To determine the effect of these properties on numerical results, a sensitivity analysis on the most important parameters affecting the fluid dynamics of a microfluidic multiphase system is required. The aim of this thesis is to compare the Volume of Fluid (VOF) and Phase Field methods implemented in OpenFOAM package, in solvers interFoam and phaseFieldFoam, respectively, in terms of their basic suitability to simulate multiphase flows. For this purpose, interFoam in foam-extend-3.1 and phaseFieldFoam in OpenFOAM-1.6-ext are used. By doing so, problematic parameters of microfluidic multiphase simulations can be identified. The case study of interest is a two dimensional motionless air bubble in stagnant water. This configuration is frequently discussed in literature as a validation ... mehr case. The simplicity of this case makes it suitable as a fast benchmark to achieve the target of this study. The starting point of the research is therefore the comparison between the numerical results obtained for this configuration with the ones reported in literature. For this benchmark problem, mismatch between discretised surface tension force and pressure gradient and/or inaccurate computation of interface curvature can lead to numerical errors, which manifest as parasitic currents. These currents strengthen when there exist discontinuities in field variables (e.g., density and pressure) at the interface between the air bubble and water. Thus, in the second step, this error and its consequences such as motion of bubble and distortion of its interface are characterised. To do so, several parameters such as spatial and temporal resolutions of the numerical solvers, initial estimation for bubble interface thickness and bubble diameter have been changed. The accuracy of the solvers is determined based on their ability to predict the pressure difference between the gas inside the bubble and the surrounding liquid and minimising the induced parasitic current and bubble motion, while maintaining an appropriate bubble interface thickness. The results show that the phaseFieldFoam is more accurate in predicting most of the characteristic properties of the considered multiphase system, compared to interFoam. Under mesh refinement, the former code converges with equal densities and viscosities but not for air-water system where the magnitude of spurious currents is independent from grid size. Furthermore, it is revealed that phaseFieldFoam requires at least six computational cells across the thickness of bubble interface to provide accurate results. It is also shown that to achieve convergence, satisfying common criteria on the group of involved non-dimensional numbers (here, Courant and Cahn numbers) is not sufficient. A further restriction originates from the computational time step. By refining the time step, the diffusion number, which is shown to influence the solution significantly, can be controlled and as a result, convergence of the results can be obtained. Moreover, the main drawback of phaseFieldFoam is illustrated to be in conserving the mass of the gas phase inside the bubble as the solver tends to shrink the bubble size. This problem, however, could be overcome by a higher mesh resolution. On the other hand, although interFoam does not suffer from problems in mass conservation, it is out of step in prediction of the characteristic parameters as well as in eliminating the parasitic current and unfavourable bubble motion. The adopted approach for tackling this problem is beneficial for numerical simulation of three-dimensional microfluidic systems with complex flow geometries, both for design and optimisation purposes.