  We develop a rapid and accurate contour method for the solution of
time-fractional PDEs. The method inverts the Laplace transform via an optimised
stable quadrature rule, suitable for infinite-dimensional operators, whose
error decreases like $\exp(-cN/\log(N))$ for $N$ quadrature points. The method
is parallisable, avoids having to resolve singularities of the solution as
$t\downarrow 0$, and avoids the large memory consumption that can be a
challenge for time-stepping methods applied to time-fractional PDEs. The ODEs
resulting from quadrature are solved using adaptive sparse spectral methods
that converge exponentially with optimal linear complexity. These solutions of
ODEs are reused for different times. We provide a complete analysis of our
approach for fractional beam equations used to model small-amplitude vibration
of viscoelastic materials with a fractional Kelvin-Voigt stress-strain
relationship. We calculate the system's energy evolution over time and the
surface deformation in cases of both constant and non-constant viscoelastic
parameters. An infinite-dimensional ``solve-then-discretise'' approach
considerably simplifies the analysis, which studies the generalisation of the
numerical range of a quasi-linearisation of a suitable operator pencil. This
allows us to build an efficient algorithm with explicit error control. The
approach can be readily adapted to other time-fractional PDEs and is not
constrained to fractional parameters in the range $0<\nu<1$.
