In order to increase the efficiency and reliability of the dynamic analysis for flexible planar linkage containing the coupling of multi-energy domains, a method based on bond graph is introduced. From the viewpoint of power conservation, the peculiar property of bond graph multiport element MTF is discussed. The procedure of modeling planar flexible multibody mechanical systems by bond graphs and its dynamic principle are described. To overcome the algebraic difficulty brought by differential causality and nonlinear junction structure, the constraint forces at joints can be considered as unknown effort sources and added to the corresponding O-junctions of system bond graph model. As a result, the automatic modeling on a computer is realized. The validity of the procedure is illustrated by a practical example.