We define a discrete-time Markov chain for abstract polymer models and show that under sufficient decay of the polymer weights, this chain mixes rapidly. We apply this Markov chain to polymer models derived from the hard-core and ferromagnetic Potts models on bounded-degree (bipartite) expander graphs. In this setting, Jenssen, Keevash and Perkins (2019) recently gave an FPTAS and an efficient sampling algorithm at sufficiently high fugacity and low temperature respectively. Their method is based on using the cluster expansion to obtain a complex zero-free region for the partition function of a polymer model, and then approximating this partition function using the polynomial interpolation method of Barvinok. Our approach via the polymer model Markov chain circumvents the zero-free analysis and the generalization to complex parameters, and leads to a sampling algorithm with a fast running time of $O(n \log n)$ for the Potts model and $O(n^2 \log n)$ for the hard-core model, in contrast to typical running times of $n^{O(\log Δ)}$ for algorithms based on Barvinok's polynomial interpolation method on graphs of maximum degree $Δ$. We finally combine our results for the hard-core and ferromagnetic Potts models with standard Markov chain comparison tools to obtain polynomial mixing time for the usual spin Glauber dynamics restricted to even and odd or `red' dominant portions of the respective state spaces.