We present an \(L^2\) norm convergence of the implicit-explicit BDF2 scheme with variable-step for the unsteady incompressible Navier-Stokes equations with an inf-sup stable FEM for the space discretization. Under a weak step-ratio constraint \(0<r_k:=\tau _k/\tau _{k-1}<4.864\) , our error estimate is mesh-robust in the sense that it completely removes the possibly unbounded quantities, such as \(\Gamma _N=\sum _{k=1}^{N-2}\max \{0,r_{k}-r_{k+2}\}\) and \(\Lambda _N=\sum _{k=1}^{N-1}(|r_{k}-1|+|r_{k+1}-1|)\) included in previous studies. In this analysis, we integrate our recent theoretical framework that employs discrete orthogonal convolution (DOC) kernels with an auxiliary Stokes problem to split the convergence analysis into two distinct parts. In the first part, we address intricate consistency error estimates for the velocity, pressure and nonlinear convection term. The resulting estimates allow us to utilize the conventional methodologies within the DOC framework to preserve spatial accuracy. In the second part, through the use of the DOC technique, we prove that the proposed variable-step BDF2 scheme is of second-order accuracy in time with respect to the \(L^2\) norm. Extensive numerical simulations coupled with an adaptive time-stepping algorithm are performed to show the accuracy and efficiency of the proposed variable-step method for the incompressible flows.