Abstarct: Viscous fingering instability is a well-known instability that occurs in the interfaces of two fluids due to the difference in viscosities. This instability could be seen in a wide variety of natural systems and industrial processes. In recent decades, its role in the process of enhanced oil recovery is perhaps the most important reason for much attention to this phenomenon. In order to access the remaining oil in the reservoir, another fluid, such as water is pumped into it, and thus residual oil can be extracted. In other hands, when a low viscous fluid (such as water) is injected into a porous medium that saturated with a high viscous fluid (such as oil), the interface of two fluids will gradually become unstable and the injected fluid penetrates with finger-like patterns into the high viscous fluid. Reaching these fingers into the main reservoirs, the oil recovery will practically stop. Therefore, a better understanding of this instability will have technological significance in order to eliminate or control it. Using polymer solution as displacing fluid is one of the useful tools for decreasing intensity of fingering instability. Typically, the polymer solutions can be classified as viscoelastic fluids. The main purpose of this study is a numerical simulation of polymer flooding and indeed, investigating fingering instability in displacements involving viscoelastic fluid. In this study, three types of the constitutive equations are used to model viscoelastic fluids: Oldroyd-B model, White-Metzner and Giesekus models. The displacements of Newtonian fluid by viscoelastic fluid through various porous media including homogeneous and isotropic medium, anisotropic medium, anisotropic velocity- dependent porous medium and heterogeneous porous media (laxyered and checkerboard systems) are simulated and the effect of different parameters of flow and media on fingering instability is studied. In linear stability analysis with the quasi-steady state approximation and shooting method, the growth rate of disturbances in different wave numbers with presented parameters is investigated. In non-linear simulation section with the accurate spectral method and taking advantages of Hartley transform, some mechanisms of growth of fingers are analyzed in different time sequences. Furthermore, an investigation of transversely average concentration, mixing length and sweep efficiency field is carried out. The results show that the elasticity of viscoelastic phase has a stabilizing effect while the shear-thinning property has a destabilizing effect. In addition, increasing permeability of media in the longitudinal direction to the transverse leads to more stable flow. On the other hand, when the ratio of transverse to longitudinal dispersion is increased, the intensity of instability decreases. When a periodic heterogeneity is considered for displacement media, the channeling regime is observed in the laxyered system and fingers interactions decrease. Increasing the number of laxyers in each direction leads to increase in sweep efficiency. In addition, the mixing length is enhanced by increasing in permeability variance in each laxyer. However, the sweep efficiency increases by increasing permeability variance due to increase breakthrough time (the time that fingers reach to end of computational domain).