A self-consistent numerical treatment for modeling fractal aggregate dynamics is presented. Fractal aggregates play an important role in a number of complex astrophysical regimes, including the early solar nebula and the interstellar medium. Aggregates can be of various forms and sizes, ranging from tiny dust particles to ice chunks in planetary rings and possibly even comets. Many observable properties, such as light scattering and polarization, may depend sensitively on the geometry and motion of such aggregates. Up to now various statistical methods have been used to model the growth and interaction of aggregates. The method presented here is unique in that a full treatment of rigid body dynamics—including rotation—is incorporated, allowing individual particle and cluster trajectories and orientations to be followed explicitly. The method involves solving Euler's equations for rigid body motion and introducing a technique for handling oblique collisions between arbitrarily shaped aggregates. Individual particles may be of any size and can have their own spin. Currently tangential impulses during impacts are assumed negligible, although equations for the treatment of tangential friction are presented. Models for the coagulation and restitution of aggregates are discussed in detail. Some of the key features required for a fragmentation model, not implemented here, are discussed briefly. Torque effects arising from self-gravity, tidal fields, or gas drag are not presently considered. Although the discussion focuses mainly on the theory behind the numerical technique, test simulations are presented to compare with an analytic solution of the coagulation equation and to illustrate the important aspects of the method.