An initial-boundary value problem for the velocity distribution of a viscoelastic flow with generalized fractional Oldroyd-B constitutive model is studied. The model contains two Riemann-Liouville fractional derivatives in time. The eigenfunction expansion of the solution is constructed. The behavior of the time-dependent components of the solution is studied and the results are used to establish convergence of the series under some conditions. Further, applying the convolutional calculus approach proposed by Dimovski (I.H. Dimovski, Convolutional Calculus, Kluwer, Dordrecht (1990)), a Duhamel-type representation of the solution is found, containing two convolution products of particular solutions and the given initial and source functions. A non-classical convolution with respect to spatial variable is used. The obtained representation is applied for numerical computation of the solution in the case of a generalized second grade fluid. Numerical results for several one-dimensional examples are given and the present technique is compared to a finite difference method in terms of efficiency, accuracy, and CPU time.