A robust algorithm for direct numerical simulation of fully three-dimensional, incompressible two-phase flow is presented. The method is introduced in the context of gas bubbles rising in viscous liquids, e.g. air bubbles rising in water. Key strengths of the simulation approach include the ability to simulate flows in an extended, wide range of Reynolds and Bond numbers and also the ability to handle large, realistic ratios of the density and viscosity of the fluids. The aptitude of the algorithm is due to the combination of several powerful components: The interface between the phases is tracked explicitly by an unstructured, adaptive, triangular mesh (front tracking), while the equations governing incompressible, Newtonian flows are solved using a finite volume scheme based on the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm. Further, the SIMPLE flow solver is integrated with PARAMESH: a block-based, adaptive mesh refinement (AMR) tool for multi-processing, allowing the solution of the governing equations in parallel on a supercomputer. Finally, the use of a non-inertial, moving reference frame attached to the rising bubble facilitates long-time simulations without requiring a larger computational domain. The methodology just described has been applied to simulate 3D gas bubbles rising in viscous liquids. In particular, air bubbles rising in water were studied, and the numerically predicted terminal rise velocities agree well with experimental measurements for bubbles with initial diameters ranging from 0.5 rum to 30 mm. Simulation results also reveal that bubbles rising in viscous liquids at high Reynolds numbers may rise on a zigzag or/and a spiral path. At last, we tried to use this model to simulate the interaction between two bubbles in liquid. The simulation results provide us with more physical insights into the complex bubble rising behavior in viscous liquids.