Based on the principle of virtual works, a multiphase smoothed particle hydrodynamics (SPH) model is further developed from the foundation of Hu X Y et al. (2006) and Grenier N et al. (2009). In the present model, the surface tension force implementation suitable for the multiphase flows with a large density ratio is applied, and this allows a good continuity at the multiphase interface. Artificial displacement correction is applied to keep the particles distributing uniformly in the whole flow field, and therefore any artificial viscous term is never needed; this is very important in the numerical simulation of viscous flows since the introduction of artificial viscosity changes the Reynolds number. Background pressure and interface sharpness force are added in the equation of state and the equation of momentum respectively to ensure the multiphase interface stability and smoothness; this is essential in the simulation of multiphase flows with large density difference at the multiphase interface. Two types of viscosity expressions suitable for multiphase flows are introduced and analyzed; the conclusion is that the formula proposed by Morris et al. (1997) and its similarly derived forms can give more accurate results. In the numerical validations, an oscillating droplet test is applied first to confirm the accuracy of the surface tension model and good results are achieved. This demonstrates that the artificial displacement and the interface sharp force will make negligible effects to the surface tension implementation. After that, two classic quantitative benchmarks of rising bubbles are simulated and the results of SPH agree well with the reference data. Moreover, in the two numerical benchmarks, the effect of the artificial displacement, the choice of the viscosity expression, and the type of the kernel function are compared and finally an optimal combination of these numerical aspects is recommended. Based on the above numerical investigations, the splitting process of an initially circular bubble is simulated and the numerical results agree well with the experimental data. In the last numerical case, the process of chasing and merging between two rising bubbles in vertical direction is simulated, based on which the mechanisms of these interesting interactions between two rising bubbles are analyzed. It is demonstrated in the present work that further improved multiphase SPH model may provide a potential method for the research of bubble dynamics.