The work presented here considers conservation and balance laws and constitutive theories for internal polar non-classical isotropic, homogeneous incompressible thermofluids presented by Surana et.al to present numerical studies and comparison with the results obtained using classical thermodynamic frame and standard constitutive theories. The internal polar continuum theories are based on the fact that if the velocity gradient tensor is a fundamental measure of deformation physics in fluids then the thermodynamic framework for such fluids must incorporate the velocity gradient tensor in its entirety. Polar decomposition of the velocity gradient tensor into stretch rates and the rotation rates shows that only the stretch rates are incorporated in the currently used thermodynamic framework and the rotation rates are completely neglected. If the velocity gradient tensor varies from a material point to the neighboring material points, then so do the rates of rotations which, when resisted by the fluid result in conjugate moment tensor. Rates of rotations and conjugate moment tensor can result in additional resistance to fluid motion and additional dissipation i.e. entropy production. Due to the fact that the internal polar non-classical continuum theory incorporates internal rotations and conjugate moment tensor, the theory is called internal polar non-classical continuum theory. The thermodynamic framework for internal polar thermofluids has been presented by Surana et.al. The constitutive theory for internal polar incompressible thermofluids has also been presented by Surana et.al. These are utilized in this work to present numerical studies for model problems. Boundary value problems consisting of fully developed flow between parallel plates, square and rectangular lid driven cavities and asymmetric sudden expansion with three different expansion ratios are used as model problems. Numerical solutions are computed using least squares finite element processes based on residual functional in which p-version hierarchical local approximations are considered in scalar product spaces that permit higher order global differentiability local approximations. Nonlinear algebraic equations resulting from the finite element formulation are solved using Newton’s linear method with line search. Numerical solutions obtained from internal polar mathematical models are compared with those obtained using classical continuum theory. Advisors/Committee Members: Surana, Karan S (advisor), TenPas, Peter W (cmtemember), Sorem, Robert M (cmtemember).