A central problem in fluid dynamics governed by the incompressible Navier-Stokes equations is verifying the global (nonlinear) stability of laminar flows against arbitrary perturbations. A well-established approach for tackling this problem is the energy method, pioneered by Reynolds and Orr. In this paper, we propose a numerical methodology in primal variables (i.e. three components of velocity and pressure) grounded in the mixed finite element (FE) analysis to study energy stability in ducts with arbitrary cross-sections. The method is sparsity-exploiting, enabling an efficient solution in domains requiring fine discretizations, such as those with small aspect ratios or with localized solution features. We present energy stability results for the fluid flows in rectangular and elliptic pipes, validated against the existing literature. Furthermore, we explore domains such as an equilateral triangle, reporting for the first time the energy stability limit of such flow, and a homotopy from a circle to a square, achieved by progressively rounding the square’s vertices. We examine the behavior of the most energy unstable mode, the energy stability limit and their variation with changes in the aspect ratio or key geometric parameter, reporting, with accuracy, any crossing in branches of eigenfunctions and whether such functions are inherently two-dimensional or three-dimensional. Full agreement between our results and existing literature reinforces the validity of our method. However, our method is more practical and general from a computational standpoint, when compared to existing methods, most notably pseudospectral methods (used in most computations), which are limited to simple geometries and provide little benefits when the solution lacks regularity. More importantly, this generality allows to compute energy stability of ducts used in practice for many engineering applications, which makes it useful in cases where a laminar flow is desired.