In contrast to our understanding of density field tracers, the modeling of direct statistics pertaining to the cosmic velocity field remains open to significant opportunities for improvement. The lack of accurate modeling for the nonlinear domain of pairwise velocities restricts our capacity to fully exploit the information encoded in this observable. We present a robust approach for modeling the mean infall velocities, v12(r ,a ), with broad applicability spanning sub-megaparsec scales and cosmologies extending beyond the standard Λ CDM paradigm. Our approach involves solving the full pair-conservation equation using accurate nonlinear power spectrum descriptions. To assess the robustness of our model, we extend it to cosmologies beyond the standard Λ CDM , in particular, the Hu-Sawicki f (R )-gravity and Dvali-Gabadadze-Porrati (DGP) modified gravity models. Remarkably, our predictions for pairwise velocities of dark matter particles at kiloparsec scales exhibit excellent agreement with N -body simulations throughout the entire dynamical range (0.1 ≲ξ ≲1000 , or r ≥0.4 h-1 Mpc ). Furthermore we show that different gravity models leave distinct signatures in the shape and dynamics of the mean pairwise velocities, providing a potent test of cosmological gravity laws.