<p>Tolerances are vital in determining the static and dynamic behavior of mechanical systems. Overly tight tolerances can be expensive or even impossible to produce, while excessively loose tolerances may introduce unwanted stresses, increasing the risk of part failure or assembly challenges. This can lead to higher scrap rates and inefficiencies, as well as noise and vibration issues that degrade product performance and durability. Despite their importance, current tolerancing methods often overlook these factors during the design stage. Existing system-level tolerance analysis is based on geometric analysis, and fails to capture more complex coupled effects. In this work, we propose a flexible multibody simulation-based approach to assess the impact of positional and orientational deviations as features of size on assembled products. In a first stage, a mesh-morphing framework is set up in order to apply geometric deviations to a nominal geometry in accordance to allocated tolerance zones. The assembly of the warped components is performed through a Flexible Natural Coordinates Formulation (FNCF). This enables an efficient analysis due to the compact load reduction and efficient joint description. Finally, this model is included in a Polynomial Chaos Expansion (PCE) to obtain the probability distribution of the performance criteria with respect to the tolerances. Case studies validate that stress distributions are highly sensitive to geometric deviation characteristics and kinematic constraints. This approach offers potential for computationally efficient dynamic simulations, sensitivity analysis, and topology optimization frameworks.</p>