The Navier-Stokes (N-S) equations for incompressible fluid flow comprise of a system of four nonlinear equations with five flow fields such as pressure P, density ρ and three velocity components u, v, and w. The system of equations is generally complex due to the fact that it is nonlinear and a mixture of the three classes of partial differential equations (PDEs) each with distinct solution methods. The N-S equations fully describe the unsteady fluid flow behaviour of laminar and turbulent types. Previous studies have shown existence of general solutions of fluid flow models but little has been done on numerical solution for velocity of flow in N-S equation of incompressible fluid flow by Crank-Nicolson implicit scheme. In practice, real fluid flows are compressible due to the inevitable variations in density caused by temperature changes and other physical factors. Numerical approximations of the general system of Navier-Stokes equations were made to develop numerical solution model for incompressible fluid flow. Adequate solutions of the latter produce numerical solutions applicable in numerical simulation of fluid flows useful in engineering and science. Non-dimensionalization of variables involved was done. Crank-Nicolson (C.N) implicit scheme was implemented to discretize partial derivatives and appropriate approximation made at the boundaries yielded a linear system of N-S equations model. The linear numerical system was then expressed in matrix form for computation of velocity field by Computational fluid dynamics (CFD) approach using MATLAB software. Numerical results for velocity field in two dimensional space, u(x,y,t) and v(x,y,t) generated in uniform 32×32 grids points of the square flow domains, 0≤x≤1.0 and 0≤y≤1.0 were presented in three dimensional figures. Results showed that the velocity in two dimensional space does not change suddenly for any change in spatial levels, x and y. Therefore, C-N implicit Scheme applied to solve the N-S equations for fluid flow is consistent.