Automated algorithms for the dynamic analysis and simulation of constrained multibody systems assume that the constraint equations are linearly independent. During the motion, when the system is at a singular configuration, the constraint Jacobian matrix possesses less than full rank and hence it results in singularities. This occurs when the direction of a constraint coincides with the direction of the lost degree of freedom. In this paper the constraint equations for deformable bodies are modified for use in the neighborhood of the singular configuration to yield the system inertia matrix which is nonsingular and also to take the actual generalized constraint forces into account. The procedures developed are applicable to both the augmented approach and the coordinate reduction methods. For the modeling of the constrained flexible multibody systems, a general recursive formulation is developed using Kane's equations, finite element method and modal analysis techniques. The system may contain revolute, prismatic, spherical or other types of joints, as well as geometrical nonlinearities; the rotary inertia is also automatically included. Simulation of a two-link flexible manipulator is presented at a singular configuration to demonstrate the utility of the method.