It has been challenging to accurately predict the unique characteristics of magnetorheological (MR) dampers, due to their inherent non-linear nature. Multidimensional flow simulation has received increasing attentions because it serves as a general methodology for modeling arbitrary MR devices. However, the compressibility of MR fluid which greatly affects the hysteretic behavior of an MR damper is neglected in previous multidimensional flow studies. This paper presents a two-dimensional (2D) axisymmetric flow of the compressible Herschel-Bulkley fluid in MR dampers. We simulated the fully coupled inertial-viscous-frictional-elastic transients in MR dampers under low-, medium-, and high frequency excitations. An arbitrary Lagrangian-Eulerian kinematical description is adopted, with the piston movements represented by the moving boundaries. The viscoplasticity and compressibility of MR fluid are, respectively, modeled by the modified Herschel-Bulkley model and the Tait equation. The streamline-upwind Petrov-Galerkin finite element method is used to solve the model equations including the conservation laws and mesh motion equation. We tested the performances of an MR damper under different electric currents and different frequency displacement excitations, and the model predictions agree well with the experimental data. Results showed that the coupled transients of an MR damper are frequency dependent. The weak compressibility of MR fluid, which mainly happens in the chamber rather than in the working gap, is crucial for accurate predictions. A damper's transition from the pre-yield to the post-yield is essentially a step-response of a second order mass-spring-viscous system, and we give such step-response a detailed explanation in terms of mass flow rate.