In this work we study a dynamical system with a complicated nonlinearity, which describes oscillation of a turbine rotor, and give an algorithm to compute Melnikov functions for analysis of its chaotic behavior. We first derive the rotor model whose nonlinear term brings difficulties to investigating the distribution and qualitative properties of its equilibria. This nonlinear model provides a typical example of a system for which the homoclinic and heteroclinic orbits cannot be analytically determined. In order to apply Melnikov's method to make clear the underlying conditions for chaotic motion, we present a generic algorithm that provides a systematic procedure to compute Melnikov functions numerically. Substantial analysis is done so that the numerical approximation precision at each phase of the computation can be guaranteed. Using the algorithm developed in this paper, it is straightforward to obtain a sufficient condition for chaotic motion under damping and periodic external excitation, whenever the rotor parameters are given.