The Navier-Stokes equations can be expressed in terms of the primary variables (e.g., velocities and pressure), secondary variables (velocity gradients, vorticity, stream function, stresses, etc.), or a combination of the two. The Least-Squares formulations of the original partial differential equations (PDE's) in terms of primary variables require C1 continuity of the finite element spaces across inter-element boundaries. This higherorder continuity requirement for PDE's in primary variables is a setback to Least-Squares formulation when compared to the weak form Galerkin formulation. To overcome this requirement, the PDE or PDE's are first transformed into an equivalent lower order system by introducing additional independent variables, sometimes termed auxiliary variables, and then formulating the Least-Squares model based on the equivalent lower order system. These additional variables can be selected to represent physically meaningful variables, e.g., fluxes, stresses or rotations, and can be directly approximated in the model. Using these auxiliary variables, different alternative Least-Squares finite element models are developed and investigated. In this research, the vorticity and stress based alternative Least-Squares finite element formulations of Navier-Stokes equations are developed and are verified with the benchmark problems. The Least-Squares formulations are developed for both the Newtonian and non-Newtonian fluids (based on the Power-Law model) and the effects of linearization before and after minimization are investigated using the benchmark problems. For the non-Newtonian fluids both the shear thinning and shear thickening fluids have been studied by varying the Power-Law index from 0.25 to 1.5. Also, the traditional weak form based penalty method is formulated for the non-Newtonian case and the results are compared with the Least-Squares formulation. The results matched with the benchmark problems for Newtonian and non- Newtonian fluids, irrespective of the formulation. There was no effect of linerization in the case of Newtonian fluids. However for non-Newtonian fluids, there was some tangible effect of linearization on the accuracy of the solution. The effect was more pronounced for lower power-law indices compared to higher power-law indices. And there seemed to have some kind of locking that caused the matrices to be ill-conditioned especially for lower values of power-law indices.