Railway ballast simulation algorithms are presented in the paper. The ballast is a system of planar rigid bodies with both convex and non-convex shapes. Contacts of interacting bodies are computed as force elements, which consist of a viscous-elastic normal part and a dry friction tangential part with sliding and sticking modes. Collision detection consists of two levels: neighbor and far ones. The far collision level detects contacts of polygon hulls by the linked linear list method. The neighbor collision level is an approach for detecting polygon penetrations by the sensitivity cell method. Simplified Jacobian matrices are used to accelerate the integration of stiff equations of motion. The ballast model can include up to some thousands of bodies and allows simulating processes of the ballast laying, compaction and so on. Some simulation results are presented.