The objective of this paper is to identify a numerical method to simulate motion of a packed or fluidized bed of fuel particles in combustion chambers, such as a grate furnace and a rotary kiln. Therefore, the various numerical methods applied in the areas of granular matter and molecular dynamics were reviewed extensively. As a result, a time driven approach was found to be suited for the numerical simulation of particle motion in combustion chambers. Furthermore, this method can also be employed to moving boundaries which are required for the present application e.g. travelling grate. The method works in a Lagrangian frame of reference, which uses the position and orientation of particles as independent variables. These are obtained by time integration of the three-dimensional dynamics equations derived from the classical Newtonian approach for each particle. This includes the keeping track of all forces and momentums acting on each particle at every time step. Viscoelastic contact forces include normal and tangential components with viscoelastic models for energy dissipation and friction. The particle shapes are approximated by spheres and ellipsoids with a varying size and ratio of the semi-axis accounting for the variety of particle geometries in a combustion chamber. For these shapes the overlap of particles during contact is expressed by a polynomial of 4th order in the two-dimensional case and a polynomial of 6th order in the three-dimensional case. A new algorithm to detect two-dimensional elliptical particle contact with sufficient accuracy was developed. It is based on a sequence of coordinate transformations and has demonstrated its reliability in numerous applications. Finally, the method was applied to simulate the motion of spherical and elliptical particles in a rectangular enclosure, on a travelling grate, and in a rotary kiln.