As an illustration for the efficiency and the robustness of our new algorithm, we applied it to the solution of a classical benchmark problem: the calculation of the flow in a cavity.

The cavity domain is the unit square and the flow is driven by a tangential velocity field applied to the top square boundary (y=1) in the absence of other body forces. Since all forces are independent of time, the flow in this case is steady and is modeled by the steady state Navier-Stokes equations with appropriate boundary conditions. The boundary conditions for the velocity are zero everywhere except on the segement y=1, where the velocity vector is (1,0).

The finite element discretization is defined in the following way. We partition the domain into (2n x 2n) square shaped elements, where n is a positive integer and define h=1/2n. Each of the square elements is further partitioned into two triangles by connecting the lower left corner to the upper right corner. The velocity space (for each component) is defined to be the space of functions that vanish on the boundary and are continuous and piecewise linear with respect to the triangulation just defined. The pressure space consists of special piecewise constant functions with respect to the square cells with zero mean value. Please refer to the paper for more details about the finite element discretization.

The iterative method is provided by an outer implicit Picard iteration combined with our algorithm for solving the corresponding nonsymmetric saddle point problem at each Picard step.

The computational results for h=1/64 and a set of different Reynolds numbers (Re=1, 10, 100, 1000) are shown below. In particular, the streamlines of the velocity field computed using this algorithm are shown. The effect of the Reynolds number on the flow pattern is clearly seen there. The flow for low Reynolds numbers (see the cases of Re=1 and Re=10) has only one vortex center, located above the center of the domain (its location moves to the right as Re increases). As Re increases further, a second vortex center appears near the lower right corner (see the case of Re=100) and, for even larger Reynolds numbers, a third vortex center develops near the lower left corner of the domain (see the case of Re=1000).

In conclusion, the implicit Picard iteration combined with our new algorithm for solving nonsymmetric saddle point problems result in a simple, robust and efficient method for solving Navier-Stokes equations for a wide range of Reynolds numbers. For each nonlinear iteration it requires the solution of a nonsymmetric saddle point problem which can be solved effectively. An advantage of this method is that it solves the discrete system without the need for additional stabilization terms in contrast to the class of penalty algorithms. The typical penalty methods add stabilization terms to the system. The bigger these terms are, the easier to solve the corresponding system is. However, these stabilization terms change the discrete equations that one solves. In particular, their presence effectively reduces the Reynolds number for the corresponding flow causing different flow behavior to be computed. On the other hand, such a problem does not exisits with the implicit method because it does not need any additional stabilization terms. The convergence of the linear solve at each Picard iteration is guaranteed only by the appropriate scaling of the preconditioners involved and by the choice of the iteration parameters (please refer to the paper for details).