This paper presents a second-order projection method for solving the time-dependent, incompressible Navier–Stokes equations. The method is based on a two-step process: first, solving diffusion-convection equations to predict intermediate velocities, and then projecting these velocities onto the space of divergence-free vector fields. The method incorporates a specialized second-order Godunov method for handling nonlinear convective terms, which provides robustness at high Reynolds numbers. The projection is approximated using a Galerkin procedure with a local basis for discretely divergence-free vector fields. The method is validated through numerical results showing second-order convergence for smooth flows and is applied to doubly periodic shear layers to assess its performance on more complex problems. The method is second-order accurate in both space and time when the time step is proportional to the spatial grid size. The algorithm is designed to handle homogeneous Dirichlet boundary conditions and assumes no external forces. The temporal discretization is based on an iterative procedure that converges to a modified Crank–Nicholson scheme, ensuring second-order accuracy. The method is tested on various flow problems, including Stokes flow, Reynolds number 100, and incompressible Euler equations, demonstrating its effectiveness and robustness. The numerical results show second-order convergence for velocity and first-order convergence for pressure, with the method being able to handle discontinuities without introducing spurious oscillations or instabilities. The method is also shown to be effective in the presence of boundaries, with the boundary layer effects being minimized through appropriate numerical techniques. The overall algorithm is efficient and well-suited for solving incompressible flow problems.This paper presents a second-order projection method for solving the time-dependent, incompressible Navier–Stokes equations. The method is based on a two-step process: first, solving diffusion-convection equations to predict intermediate velocities, and then projecting these velocities onto the space of divergence-free vector fields. The method incorporates a specialized second-order Godunov method for handling nonlinear convective terms, which provides robustness at high Reynolds numbers. The projection is approximated using a Galerkin procedure with a local basis for discretely divergence-free vector fields. The method is validated through numerical results showing second-order convergence for smooth flows and is applied to doubly periodic shear layers to assess its performance on more complex problems. The method is second-order accurate in both space and time when the time step is proportional to the spatial grid size. The algorithm is designed to handle homogeneous Dirichlet boundary conditions and assumes no external forces. The temporal discretization is based on an iterative procedure that converges to a modified Crank–Nicholson scheme, ensuring second-order accuracy. The method is tested on various flow problems, including Stokes flow, Reynolds number 100, and incompressible Euler equations, demonstrating its effectiveness and robustness. The numerical results show second-order convergence for velocity and first-order convergence for pressure, with the method being able to handle discontinuities without introducing spurious oscillations or instabilities. The method is also shown to be effective in the presence of boundaries, with the boundary layer effects being minimized through appropriate numerical techniques. The overall algorithm is efficient and well-suited for solving incompressible flow problems.