Numerical simulations by the volume of fluid model (3D-VOF) were performed for the dynamical behavior of bubbles with equilateral triangle arrangement continuously rising in shear-thinning fluids. The influences of bubble diameter, bubble formation frequency, initial bubble distance and liquid rheological property on the interaction between bubbles and drag coefficient were studied systematically. With the increase of bubble diameter, bubble formation frequency and the decrease of initial bubble distance, both the horizontal attractive effect and vertical wake effect of bubbles strengthened, resulting in the decrease of drag coefficient. The intensified liquid shear-thinning property for the lower flow index could also enhance these two effects, and accordingly lead to a decrease of bubble drag coefficient. To characterize the interaction within multi-bubble systems, the bubble interaction coefficient Ψ was introduced, and an empirical correlation of Ψ was proposed. Based on the Magnaudet model, a new bubble drag coefficient relationship for the multi-bubble system was proposed by Ψ with an average deviation of 5%.