Currently existing computational fluid dynamics-discrete element method (CFD-DEM) solvers suffer from computationally expensive coupling between the CFD and DEM as it requires calculating at each fluid time-step the void fraction and the solid-fluid forces such as drag, lift, buoyancy, and undisturbed flow forces. We develop a unified finite element CFD-DEM solver which integrates the CFD and DEM solvers into a single software resulting in faster and cheaper coupling. Our fluid formulation is stabilized using tailored techniques to prevent oscillations in regions of sharp gradients, to enhance the robustness of the formulation and local mass conservation, and to relax the Ladyzhenskaya-Babuska-Brezzi inf-sup condition. The developed solver supports high order finite elements resulting in better accuracy with larger cell sizes. Moreover, our solver supports dynamic load balance parallelization for both the particles and the fluid. This evens the distribution of workloads among processors, resulting in better efficiency and resource exploitation. Additionally, we develop a new spatially and temporally continuous analytical void fraction scheme called Quadrature-Centered Method (QCM). This scheme results in less computational time, better accuracy and convergence, and enhanced mass conservation. It also enables the use of very small CFD time-steps thus achieving better temporal accuracy and the use of mesh sizes smaller than those commonly used in CFD-DEM (< 3 times the particle diameter.) We validate our solver through several cases among which we will discuss a spouted bed test case where particle velocities matched those of the experiments, a particle sedimentation test case where we study the effect of the void fraction scheme choice on the particle’s terminal velocity, and a particle Rayleigh-Taylor Instability where the particles constituted the heavy phase and where we study the evolution of the mixing layer with time.