Based on the variational constraint approach, the variational form of Reynolds equation in hydrodynamic lubrication is revised continuously to satisfy certain constraints in the cavitation zone of oil film field. According to the physical characteristic of oil film, an eight-node isoparametric finite element method is used to convert the revised variational form of Reynolds equation to a discrete form of finite dimensional algebraic variational equation. By this approach, a perturbance equation can be obtained directly on the finite element equation. Consequently, nonlinear oil film forces and their Jacobian matrices are calculated simultaneously, and compatible accuracy is obtained without increasing the computational costs. A method, which is a combination of predictor corrector mechanism and Newton Raphson method, is presented to calculate equilibrium position and critical speed corresponding to Hopf bifurcation point of bearing-rotor system, as by-product dynamic coefficients of bearing are obtained. The timescale, i.e., the unknown whirling period of Hopf bifurcation solution of bearing-rotor system is drawn into the iterative process using Poincar?Newton Floquet method. The stability of the Hopf bifurcation solution can be detected when estimating Hopf bifurcation solution and its periods. The nonlinear unbalanced periodic responses of the system are obtained by using PNF method and path-following technique. The local stability and bifurcation behaviors of periodic motions are analyzed by Floquet theory. Chaotic motions are analyzed by Lyapunov exponents. The numerical results revealed the rich and complex nonlinear behavior of the system, such as periodic, quasiperiodic, jumped solution, chaos, and coexistence of multisolution, and so on.