In this paper we discuss the model of fractional oscillator where the inertial and restoring force terms maintains their usual expression but the damping term involves a fractional derivative of Caputo type, the so called fractional Kelvin-Voigt oscillator. The transient solution of this model is given in terms of the so called bivariate Mittag-Leffler function, while the steady-state solution in response to a sinusoidal force involves a 4-variate Mittag-Leffler function. We give numerical examples comparing the solutions for different values of the order $\alpha$ of the fractional derivative ( $0<\alpha \le 1$), and compare them with the usual $\alpha =1$ solutions in the underdamped, overdamped and critically damped situations.