By "the simplest" I mean the simplest to learn and implement from scratch. I hope my question can more or less be answered.
In two dimensions, the velocity-vorticity formulation is the simplest to implement because the variables are collocated, but boundary conditions can be complicated and it's a less direct statement of the problem. For primitive variable formulations, the staggered grid finite difference method of Harlow and Welch (1965) is a great place to start.
You can find a fully documented implementation of a very simple, yet quite efficient, solution method (Chorin's splitting method) here.
For a selection of other popular methods, take a look at Chapter 21 in this book.
Disclaimer: I'm the (co)author of both the demo program and the book. The book can be downloaded for free.
Simplest is always going to be relative to your particular interests and needs. I agree with Anders that, for incompressible flow on domains with simple geometry, you'd be hard-pressed to beat the projection method (i.e., Chorin's splitting method) if you are prioritizing both ease of use and accuracy.
To go into a little more detail, the method is question is introduced in . The more modern, second-order, approximate projection method is explained well in . The motivation is that solving the full incompressible Navier-Stokes equations requires solving for the velocity field and the pressure simultaneously, and the resulting linear system is rather ill-conditioned. The projection method eliminates this problem by splitting each time step into a velocity solve, using the pressure from the previous timestep, followed by a pressure update, which essentially enforces that the velocity field remain incompressible.
To implement this, you'll need a few other components, but all can be learned and programmed quite easily.
For the pressure solve, assuming you're interested in systems with constant density, you'll need to solve Poisson's equation. There are, of course, dozens of algorithms to approach this problem, but by far the easiest to implement -- if perhaps not to fully understand -- is the conjugate gradient (CG) algorithm. One of the best explanations of CG I've read was written by Jonathan Shewchuk and can be found here. You certainly need not read the entire paper, however, to be able to simply implement the algorithm.
You'll need another algorithm to handle the advection term in Navier-Stokes. In several dimensions, programming robust implementations of the most flexible methods, e.g. Godunov, can be quite challenging. However, provided you're interested in flows with relatively modest Reynolds number (i.e., with non-negligible viscosity), one of the essentially non-oscillatory (ENO) methods fits the bill nicely in terms of ease of implementation. There is an excellent overview of both the theory and implementation in .
You'll need to handle the viscous term using an implicit method, typically Crank-Nicolson. This is explained in detail in the projection method papers, and you can easily use CG for the matrix solve provided viscosity is constant.
 A. J. Chorin, Numerical solution of the Navier-Stokes equations, J. Math. Comput., 22 (1968), pp. 745-762
 A. Almgren, J.B. Bell, and W. Szymczak, A numerical method for the incompressible Navier-Stokes equations based on an approximate projection, SIAM J. Sci. Comput. 17 (1996), pp. 258-369.
 S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. Springer-Verlag New York,. Applied Mathematical Sciences, 153, 2002
Computer graphics and games has seen a huge explosion of interest in fluid simulation in recent years. Here is a great paper from Jos Stam which discusses implentation of a solver for realtime applications. It comes with very easy to understand source code. I don't know how accurate it is, but it might be what you're looking for.
Another really nice and simple method is using cellular automatons for discretization. There are plenty of such models, including LBA, FHP and much more. These are really nice since they can provide a realtime simulation on modern computers and can be also nicely paralellized and ran on GPUs. They also have some disadvantages and the results are strongly dependent on the shape of lattice applied. Square lattice is insufficient because it lacks rotational freedom and i.e. von kaarman vortices will be square shaped which is not nice :)