In order to increase the efficiency and reliability of dynamic and static force analysis for spatial multibody systems containing the coupling of multi-energy domains, a method based on gyrobondgraph is introduced. The procedure of modeling spacial multibody systems by bond graphs and its dynamic principle are described. The unified formulae of driving moment and constraint forces at joints are derived based on gyrobondgraph, they are easily generated on a computer in a complete form. As a result, the very difficult algebraic problem caused by differential causality and nonlinear junction structure can be overcome, and the automatic dynamic analysis of multibody systems on a computer is realized. By a practical example, the validity of this procedure is illustrated.