This paper presents a stable penalty-based model for simulating frictional contacts between many complex objects. The major advantage of our model is that it solves the problems in implementing Coulomb's friction model for computer simulation: iterative computation and slip velocity threshold. We also introduce a robust method for computing the normal vector and penetration depth at each contact point of a pair of interpenetrating polygonal objects. We demonstrate the validity and usability of the model by comparing the simulation results with closed-form solutions of Coulomb's friction model and conventional friction model, as well as performing dynamics simulation of highly complex scenes with tens of objects composed of thousands of polygons