Magnetorheological (MR) fluids are composed of a solution of micron-sized magnetizable particles dispersed in a carrier fluid, typically oil. Thus, the MR fluids are classified as smart fluids by showing varying apparent viscosity under a precisely controllable external magnetic field. This feature of the MR fluids gave the possibility to use them in the MR dampers to control variational damping forces. In this work, an MR damper on relatively small scales is considered to investigate the MR fluid flow behavior through the piston annuli where the magnetic field exists. This work mainly focuses on the numerical solutions of the flow variables in a magnetically excited non-Newtonian flow medium by comparing analytical results for various Reynolds numbers. The Herschel-Bulkley (HB) viscous model is used for the non-Newtonian characteristic of the MR fluid. For the numerical modeling of the HB viscous model, the regularized approaches are used to avoid numerical errors. However, on the contrary to the actual HB model, the regularized models can be incapable to give true apparent viscous values at very low shear rates depending on the regularization parameters. Thus, the present study aims to give a better understanding of choosing the optimal regularization parameters for the studied flow conditions. The second part of the study discusses the numerical discretization schemes aiming to present their performances while changing the Reynolds number between 0.002-1. In this manner, a CFD solver has been developed for two-dimensional geometries which are meshed with structured grids by using six different discretization schemes including two of the most known Total Variation Diminishing (TVD) schemes. For simplicity of the problem, a two-dimensional between parallel plates geometry and a constant magnitude of the magnetic field intensity along the piston annuli is assumed. The study also summarizes the CFD technique by evaluating the physical meaning of the flow field variables.