As a universal quantum mechanical approach to the dynamical many-body problem, the time-dependent density functional theory (TDDFT) might be inadequate to describe crucial observables that rely on two-body evolution behavior, like the double-excitation probability and two-body dynamic correlation. One promising remedy is to utilize the time-dependent 2-reduced density matrix (2-RDM) that directly represents two-body observables in an N-particle system. However, the usage of 2-RDM is prohibitive due to the augmented dimensionality. This work addresses the high-dimensional numerical challenges by using an equivalent Wigner phase-space formulation of 2-RDM and seeking efficient spectral approximations to nonlocal quantum potentials. For spatial periodic case, a pseudo-difference operator approach is derived for both the Hartree-exchange-correlation term and two-body collision operator, while the discretization via the Chebyshev spectral element method is provided for non-periodic case. A massively parallel numerical scheme, which integrates these spectral approximations with a distributed characteristics method, allows us to make the first attempt for real simulations of 2-RDM dynamics. Numerical experiments demonstrate the two-body correction to the quantum kinetic theory, and show the increase in the system`s entropy induced by the two-body interaction. This may pave the way for an accurate description of 2-RDM dynamics and a practical application of many-body TDDFT.