The divergence-free finite elements introduced in this paper are derived from Hermite functions, which interpolate stream functions. Velocity bases are derived from the curl of the Hermite functions. These velocity basis functions constitute a solenoidal function space, and the gradient of the Hermite functions constitute an irrotational function space. The incompressible Navier-Stokes equation is orthogonally decomposed into its solenoidal and irrotational parts, and the decoupled Navier-Stokes equations are then projected onto their corresponding spaces to form appropriate variational formulations. The degrees of the Hermite functions we introduce in this paper are bi-cubis, quartic, and quintic. To verify the accuracy and convergence of the present method, three well-known benchmark problems are chosen. These are lid-driven cavity flow, flow over a backward facing step, and buoyancy-driven flow within a square enclosure. The numerical results show good agreement with the previously published results in all cases.