This paper introduces a new method for proving global stability of fluid flows through the construction of Lyapunov functionals. For finite dimensional approximations of fluid systems, we show how one can exploit recently developed optimization methods based on sum-of-squares decomposition to construct a Lyapunov function. We then show how these methods can be extended to full infinite dimensional Navier-Stokes systems using robust optimization techniques. Crucially, this extension requires only the solution of infinitedimensional linear eigenvalue problems and finite-dimensional sum-of-squares optimization problems. The resulting method is guaranteed to provide results at least as good as classical energy methods in determining a lower-bound on the maximum Reynolds number for which a flow is globally stable.