  We present a new development of the force-coupling method (FCM) to address
the accurate simulation of a large number of interacting micro-swimmers. Our
approach is based on the squirmer model, which we adapt to the FCM framework,
resulting in a method that is suitable for simulating semi-dilute squirmer
suspensions. Other effects, such as steric interactions, are considered with
our model. We test our method by comparing the velocity field around a single
squirmer and the pairwise interactions between two squirmers with exact
solutions to the Stokes equations and results given by other numerical methods.
We also illustrate our method's ability to describe spheroidal swimmer shapes
and biologically-relevant time-dependent swimming gaits. We detail the
numerical algorithm used to compute the hydrodynamic coupling between a large
collection ($10^4-10 ^5$) of micro-swimmers. Using this methodology, we
investigate the emergence of polar order in a suspension of squirmers and show
that for large domains, both the steady-state polar order parameter and the
growth rate of instability are independent of system size. These results
demonstrate the effectiveness of our approach to achieve near continuum-level
results, allowing for better comparison with experimental measurements while
complementing and informing continuum models.
