Simulating incompressible fluid motion remains a central topic in physics-based simulation. The goal is to construct discrete fluid motions that faithfully reproduce the dynamical phenomena of their continuous counterparts. This has motivated a long line of work into \emph{structure-preserving} discretizations of the incompressible Euler equations, which aim to discretize the geometric interpretation of incompressible Euler solutions as geodesics on an infinite-dimensional manifold of volume-preserving diffeomorphisms. The Koopman operator of such diffeomorphisms is orthogonal. Motivated by this, a line of work beginning with \citet{Pavlov:2011:SPD} seeks to spatially discretize flow maps as orthogonal matrices, and model an ideal fluid by a least action principle on the orthgonal group. To make this work, one picks a finite-dimensional space of velocity vector fields, associated to skew-symmetric matrices. The resulting spatially discrete fluid model is non-holonomically constrained. Prior work has invoked the Lagrange-d'Alembert (LdA) variational principle for non-holonomically constrained systems [Pavlov et. al 2011, Gawlik et al. 2011, Natale & Cotter 2018]; this principle effectively imposes the constraint by applying friction-like force [Abanov & Khesin, 20215] on the equations of motion, which are no longer symplectic. The primary goal of this work establishes the \emph{vakonomic} variational principle as an alternative to LdA for incompressible fluid simulations. The vakonomic principle allows only variations which obey the constraints; based on this, we obtain a new discrete equation that is a Hamiltonian flow under the Lie-Poisson structure for so(n)^*. This equation can also be understood as a \emph{sub-Riemannian geodesic}, mirroring the classical Euler-Arnold equations. This has the enormous positive consequence of maintaining symplecticity, so that Casimirs are conserved to machine precision. We also illustrate a practical means of computing this Hamiltonian flow using Clebsch variables. Finally, our system can also be interpreted as an isospectral Lax equation, much like Zeitlin's Matrix Hydrodynamics model [Modin & Viviani 2025]. Altogether, this approach is a flexible framework yielding an intriguing Hamiltonian discretization of the Euler equations.