<p>This study presents a two-dimensional fluid model for simulating DC corona discharge in a point-to-plane configuration, implemented using the Least Squares Finite Element Method (LSFEM). The model accounts for electrons, positive ions, and negative ions, with transport described by the drift–diffusion approximation and coupled self-consistently to the electric field through Poisson’s equation. Field-dependent ionization, attachment, and recombination processes are incorporated to capture the essential physics of negative corona discharges. The LSFEM framework enables stable and accurate resolution of steep gradients and complex boundary conditions, while maintaining computational efficiency over long simulation times. Results reveal the temporal and spatial evolution of charged species, the redistribution of the electric field under space-charge effects, and the transition toward a quasi-steady state dominated by negative ion accumulation. Comparisons with previous numerical and experimental studies demonstrate that the proposed approach reproduces key features of corona behavior, offering a robust and efficient tool for both physical understanding and engineering design of corona-based systems.</p>