We present an automated framework for the calculation of beam functions that describe collinear initial-state radiation at hadron colliders at next-to-next-to leading order (NNLO) in perturbation theory. By exploiting the infrared behaviour of the collinear matrix elements, we factorise the phase-space singularities with suitable observable-independent parametrisations. Our numerical approach applies to a large class of collider observables, and as a check of its validity, we compute the quark beam functions for transverse-momentum resummation and N-jettiness, which are known analytically at this order, finding excellent agreement.