Understanding how variation impacts a multibody dynamic (MBD) system's response is important to ensure the robustness of a system. However, how the variation propagates into the MBD system is complicated because MBD systems are typically governed by a system of large differential algebraic equations. This paper presents a novel process, variational work, along with the polynomial chaos multibody dynamics (PCMBoD) automation process for utilizing polynomial chaos theory (PCT) in the analysis of uncertainties in an MBD system. Variational work allows the complexity of the traditional PCT approach to be reduced. With variational work and the constrained Lagrangian formulation, the equations of motion of an MBD PCT system can be constructed using the PCMBoD automated process. To demonstrate the PCMBoD process, two examples, a mass-spring-damper and a two link slider–crank mechanism, are shown.