This study proposes a robust mathematical and computational framework for investigating the dynamic response of multidirectional functionally graded viscoelastic beams (MDFGVBs) on elastic foundations under moving loads. To capture the time-dependent behavior of viscoelastic materials, the Kelvin–Voigt constitutive model is adopted to capture both stiffness and damping effects. The beam's material properties vary continuously and independently along the axial, transverse, and thickness directions, governed by a generalized spatial gradation function, reflecting realistic material heterogeneity. The governing equations of motion are systematically derived using both Euler–Bernoulli and Timoshenko beam theories, incorporating bending, shear deformation, and rotary inertia effects. The Ritz method is employed for spatial discretization under various classical boundary conditions, while transient dynamic responses are computed using the unconditionally stable Newmark-β scheme. Model accuracy and numerical stability are validated through rigorous comparisons with established benchmark solutions. A comprehensive parametric study is conducted to explore the influence of foundation stiffness, viscoelastic damping coefficients, and material gradation indices on both free and forced vibrations. The results reveal significant coupled effects between the foundation parameters and the spatially varying viscoelastic and elastic properties. These findings provide practical guidance for designing functionally graded beam components, for instance in aerospace wing structures and precision robotic arms, where controlling vibration amplitudes and dynamic deflections is critical for structural integrity and operational accuracy. The study also aids in selecting appropriate material gradation profiles and foundation stiffness to achieve targeted dynamic performance.