A gas-liquid two phase flow is complicated and it has not been understood well thus far, in spite of extensive investigation. Numerical simulation is a potential approach to understand this phenomenon. Although a number of studies have been conducted to understand the behavior of bubbles on the basis of computational fluid dynamics (CFD), it is difficult to completely simulate a complicated three-phase flow, including coalescence and breakup of bubbles. Although the two-fluid model based on the semi-empirical model can well estimate the actual behavior of the system in which the equations are derived, the estimation over the applicable region of equations does not always agree with the actual result. Since the 1960s, various procedures have been proposed to directly track the free surface between two phases, for example, the adaptive mesh method and the particle method. Although each of these methods has certain advantages and disadvantages, the volume of fluid (VOF) method is the most acceptable method for capturing the free surface accurately and clearly. However, a concern related to this method is the maintenance of a constant volume of the fluid. In this study, a simulation code using the VOF method is developed in order to estimate the behavior of bubbles in a vertical pipe. Further, an offset of the volume fraction is introduced to stably calculate and minimize the volume fluctuation. The effect of the surface tension is also built into the program in order to estimate the behavior of the bubbles rising through the liquid medium. The simulations of the collapsed water column and a single rising bubble are conducted with the proposed simulation code. Consequently, we confirm that these results fairly agree with the experimental ones.