In this paper, we proposed a fully discrete projection method with modular grad-div stabilization for solving the time-dependent inductionless magnetohydrodynamic equations. The method incorporates a minimally intrusive module into the classical projection method, serving as a post-processing step, thereby enhancing solution accuracy and improving mass conservation. Concurrently, a decoupled strategy is employed to separate the magnetic and fluid field functions from the original system. Therefore, at each time step, we only need to solve several linear sub-systems for which the numerical solutions can be obtained efficiently. For spatial discretization, the current density and electric potential are discretized by H ( div , Ω ) × L 2 ( Ω ) -conforming finite element pair, which ensures that the discrete current density is exactly divergence-free. Therefore, the designed numerical scheme maintains the features of linearization, decoupling, unconditional energy stability, charge conservation, and improved mass conservation. The unconditional energy stability and convergence of the algorithm are analyzed and derived. Numerical results are presented to verify that the algorithm exhibit robustness with respect to the stabilization parameters and demonstrate the performance of the scheme, particularly with respect to its stability and accuracy.
