- Due to simplicity in implementation and data structure, elements with equal-order interpolation of velocity and pressure are very popular in finite-element-based flow simulations. Although such pairs are inf-sup unstable, various stabilization techniques exist to circumvent that and yield accurate approximations. The most popular one is the pressure-stabilized Petrov–Galerkin (PSPG) method, which consists of relaxing the incompressibility constraint with a weighted residual of the momentum equation. Yet, PSPG can perform poorly for low-order elements in diffusion-dominated flows, since first-order polynomial spaces are unable to approximate the second-order derivatives required for evaluating the viscous part of the stabilization term. Alternative techniques normally require additional projections or unconventional data structures. In this context, we present a novel technique that rewrites the second-order viscous term as a first-order boundary term, thereby allowing the completeDue to simplicity in implementation and data structure, elements with equal-order interpolation of velocity and pressure are very popular in finite-element-based flow simulations. Although such pairs are inf-sup unstable, various stabilization techniques exist to circumvent that and yield accurate approximations. The most popular one is the pressure-stabilized Petrov–Galerkin (PSPG) method, which consists of relaxing the incompressibility constraint with a weighted residual of the momentum equation. Yet, PSPG can perform poorly for low-order elements in diffusion-dominated flows, since first-order polynomial spaces are unable to approximate the second-order derivatives required for evaluating the viscous part of the stabilization term. Alternative techniques normally require additional projections or unconventional data structures. In this context, we present a novel technique that rewrites the second-order viscous term as a first-order boundary term, thereby allowing the complete computation of the residual even for lowest-order elements. Our method has a similar structure to standard residual-based formulations, but the stabilization term is computed globally instead of only in element interiors. This results in a scheme that does not relax incompressibility, thereby leading to improved approximations. The new method is simple to implement and accurate for a wide range of stabilization parameters, which is confirmed by various numerical examples.…