Since it is difficult to analytically express the Melnikov function when a dynamic system possesses multiple saddle fixed points with homoclinic and/or heteroclinic orbits, this paper investigates a vehicle model with nonlinear suspension spring and hysteretic damping element, which exhibits multiple heteroclinic orbits in the unperturbed system. First, an algorithm for Melnikov integrals is developed based on the Melnikov method. And then the amplitude threshold of road excitation at the onset of chaos is determined. By numerical simulation, the existence of chaos in the present system is verified via time history curves, phase portrait plots and $Poincar{acute{e}}$ maps. Finally, in order to further identify the chaotic motion of the nonlinear system, the maximal Lyapunov exponent is also adopted. The results indicate that the numerical method of estimating chaotic threshold is an effective one to complicated vehicle systems.