Adding damping devices to the rotor supports is a frequently used technological solution for reducing vibrations of rotating machines. To achieve their optimum performance, their damping effect must be adaptable to the current operating speed. This is offered by magnetorheological squeeze film dampers. The magnetorheological oils are liquids sensitive to magnetic induction and belong to the class of fluids with a yielding shear stress. Their response to the change of a magnetic field is not instantaneous, but it is a process called the delayed yielding. The developed mathematical model of the magnetorheological squeeze film damper is based on the assumptions of the classical theory of lubrication. The lubricant is represented by a bilinear material, the yielding shear stress of which depends on magnetic induction. The delayed yielding process is described by a convolution integral with an exponential kernel. The developed mathematical model of the damper was implemented in the computational procedures for transient analysis of rotors working at variable operating speed. The carried-out simulations showed that the delayed yielding effect could have a significant influence on performance of magnetorheological damping devices. The development of a novel mathematical model of a magnetorheological squeeze film damper, the representation of the magnetorheological oil by bilinear material, taking the delayed yielding phenomenon into consideration, increased numerical stability of the computational procedures for transient analysis of flexible rotors, and extension of knowledge on behavior of rotor systems damped by magnetorheological squeeze film dampers are the principal contributions of this paper.