Based on the continuous Galerkin method, a time-discrete variational integrator was presented for multibody dynamic systems with nonholonomic constraints. The weighted residual statements of Hamilton’s canonical equations were taken firstly. Then time-stepping schemes were outlined, and algorithmic conservations were discussed. Finally, a simplified model of the skateboard validated the accuracy and efficiency of the method presented.