Flow physics and heat transport properties of turbulent Rayleigh--Bénard convection with polymer additives are investigated via direct numerical simulations (DNS) in a cylindrical cell with no-slip top, bottom and lateral boundary conditions. Three-dimensional DNS are performed for a Weissenberg number range of 0 ≤Wi ≤ 50 with the FENE-P (finitely extensible nonlinear elastic-Peterlin) viscoelastic constitutive model at three Rayleigh numbers Ra = 10^7,10^8,5*10^8 and one Prandtl number Pr = 4.3. The solvent viscosity ratio is fixed at 0.9 and the maximum chain extensibility parameter L = 50, corresponding to the dilute polymer solutions with moderate fluid elasticity. The polymer additives are found to suppress the small-scaled fluctuation structures and enhance the plume coherency. At Ra = 10^7, a small heat enhancement is obtained at small Wi < 5 and followed by a monotonic heat reduction up to 7%; Remarkably, the DNS results have demonstrated that with increasing Wi the polymers stabilize the flow, gradually organizing it into an axisymmetric pattern. Specifically, the flow evolves from the initial Newtonian large-scale circulation, into transitional toroidal rings with bulk up/down-welling, and finally into a fountain-shaped pattern—a qualitative trend consistent with Keqing Xia's recent experiments. At higher Ra ≥ 10^8, the polymers cause monotonic heat and momentum transport reduction with increasing Wi, with a trend to approach a saturate heat and momentum reduction states, which also agrees qualitatively with the experiments. However, in our DNS, for Ra = 5*10^8, a heat transport reduction up to 12% is obtained, which is much smaller than the 44% reduction obtained in experiments. This discrepancy could be originated from the moderate L = 50 used in our simulations. To this end, we performed another set of DNS at Ra = 10^7 with the Oldroyd-B constitutive model, where L → ∞. As expected, with increasing Wi, polymers cause monotonic heat transport reduction up to 60%. The turbulent statistics in terms of velocity, temperature and elastic stresses are also discussed to provide more insights into determining the numerical rheological parameters for real experimental polymer solutions used in turbulent flows.