ABSTRACTThe simulation of polymer injection in a reservoir is of paramount importance in enhanced oil recovery. Despite decades of research, the computation of polymer flows in porous media remains a challenging task. The main difficulty lies in the necessity to take into account the effect of inaccessible pore volumes (IPV), for which standard closure laws give rise to a weakly hyperbolic or even non‐hyperbolic system. In the latter case, exponential instabilities may appear at the continuous level, which must be addressed at the discrete level so as to prevent a premature stop of the numerical simulations. The notion of IPV was introduced by engineers in order to account for the following observation: when a polymer solution is injected into an initial core saturated with water, the breakthrough of the polymer at the exit occurs before that of the water in which it is injected. It seems that due to their large size, the polymer molecules cannot insinuate themselves into all pores as well as water. Having less volume to flood, the polymer molecules see their speed increased, hence the ad hoc acceleration factor associated with the polymer. In this work, we propose a relaxation method that guarantees some practical robustness for all IPV laws. This is achieved by replacing the original system by a relaxation model which is always hyperbolic. The designed relaxation model involves two parameters which enable us not only to adjust the correct amount of numerical dissipation, but also to ensure positivity for some critical quantities such as water saturation and polymer concentration. Extensive numerical tests are performed in order to compare the relaxation scheme to the more classical upwind scheme for several IPV laws.