This article presents the novel algorithm using 2D Bernstein polynomials in the solution of fractional derivative partial differential equation of dynamic viscoelastic (VE) beam rested on nonlinear elastic foundations, under various loading scenarios. Based on Euler–Bernoulli thin beam theory and the fractional Kelvin–Voigt viscoelastic model, the nonlinear multi-fractional partial differential equation governing the VE beam is established. Firstly, the integer and fractional differential matrices of Bernstein polynomials in one dimensional are derived. Then, two-dimensional Bernstein polynomial operational matrices (2D-BPOM) for integer order and fractional order of differentiation is deduced. The 2D-BPOM is employed to discretize the governing nonlinear partial differential equation into a system of nonlinear algebraic equations, which are solved via Newton’s method. Verifications with exact solutions and numerical ones are presented to proof the solution technique and validate mathematical model. Comprehensive parametric studies are performed to examine the impact of loading conditions, foundation parameters and fractional orders on the dynamic response of VE beam. This study is limited to using a constant-order fractional derivative within the Kelvin–Voigt viscoelastic model for the VE thin beam’s governing equation neglecting a shear effect. The findings of this study may enable researchers to select an appropriate mathematical model that accurately aligns with a specific experimental model.