Squirmers are models of a class of microswimmers that self-propel in fluids without significant deformation of their body shape, such as ciliated organisms and phoretic particles. Available techniques for their simulation are based on the boundary element method and do not contemplate nonlinearities such as those arising from the inertia of the fluid or non-Newtonian rheology. This article describes a methodology to simulate squirmers that overcomes these limitations by using volumetric numerical methods, such as finite elements or finite volumes. It deals with interface conditions at the surface of the squirmer that generalize those in the published literature, which are generally restricted to the imposition of slip velocities. The actual procedures to be performed on a fluid solver to implement the proposed methodology are provided, including the treatment of metachronal surface waves. Among the several numerical examples, a two-dimensional simulation is shown of the hydrodynamic interaction of two individuals of Opalina ranarum.