A three-dimensional detailed numerical simulation for flows with bubbles was conducted in the present study. The simulation was based on the local instantaneous field equations of the gas liquid two-phase flow and on the interface tracking method of the volume of fluid method. The velocity and pressure fields were calculated using the modified SOLA method. The validity of the numerical method was confirmed by the comparison between measured and calculated terminal velocities and shapes of a single bubble in an infinite stagnant liquid for a wide range of the Morton number and the Eotvos number. Then, laminar bubbly flow in a vertical square duct was analyzed by making use of a periodic boundary condition. As a result, two typical void distributions, the wall peak and core peak, were obtained by varying the Eotvos number. This result agreed well with existing experimental data.