We introduce a novel dual exactly divergence-free and high-order virtual element method for solving the inductionless magnetohydrodynamic (MHD) equations on a general Lipschitz domain. The key advantage of the proposed method is its ability to naturally preserve the pointwise divergence-free properties of the velocity and current density in a discrete sense by utilizing specially designed virtual element spaces [17, 24]. To achieve higher accuracy (i.e., greater than second-order accuracy for \(k\ge 2\) ), it suffices to naturally extend the virtual element spaces and degrees of freedom to accommodate polynomial projections of degree k. We rigorously prove the well-posedness of the discrete scheme and establish its convergence analysis. It is worth noting that the error components partially decouple. Specifically, the error estimates for velocity and current density are asymptotically robust with respect to pressure and electric potential, while the error estimates for pressure and electric potential are mutually asymptotically robust. We finally present a series of numerical simulations to validate the theoretical predictions and demonstrate the effectiveness of the proposed method.
