Based on the Eulerian finite element method (EFEM), an axially symmetric numerical model was established for 2-bubble underwater pulsation. The accuracy of the model and the convergence of the mesh were fully verified by comparison with the unified bubble equation and experimental results. The calculation results show that, the unified bubble equation is more accurate than other bubble theories in predicting the bubble dynamic behavior and the pressure load in the flow field. Combined with the EFEM and the unified bubble equation, the effects of buoyancy parameter δ and strength parameter ε on the coupling law of double bubbles were studied. For buoyancy parameter δ≤0.15, the upper bubble will produce a vertical downward jet under the action of the lower bubble, and the lower bubble boundary is like the solid wall boundary. For δ increasing to 0.2, the influence of the lower bubble on the upper bubble weakens, the buoyancy effect becomes more prominent and the jet direction of the upper bubble is vertical upward. Strength parameter ε has no obvious effect on the coupling between bubbles, but its effect on the bubble jet velocity decreases significantly for ε≥150.