Antiferromagnetic spintronic devices have the potential to greatly outperform conventional ferromagnetic devices due to their ultrafast dynamics and high data density. A challenge in designing these devices is the control and detection of the orientation of the anti-ferromagnet. One of the most promising ways to achieve this is through the exchange bias effect. This is of particular importance in large scale multigranular devices. Previously, due to the large system sizes, only micromagnetic simulations of exchange have been possible, with an assumed a distribution of antiferromagnetic anisotropy directions and grain size. Here, we use an atomistic model where the distribution of antiferromagnetic anisotropy directions occurs naturally and where the exchange bias occurs due to the intrinsic disorder in the antiferromagnet. We perform large scale simulations of exchange bias, generating realistic values of exchange bias. We find a strong temperature dependence of the exchange bias in agreement with experimental observations, approaching zero at the blocking temperature of the antiferromagnet. We find that the experimentally observed increase in the coercivity at the blocking temperature occurs due to the superparamagnetic flipping of the antiferromagnet during the hysteresis loop cycle. We find a large discrepancy between the exchange bias predicted from a geometric model of the antiferromagnetic interface indicating the importance of grain edge effects in multigranular exchange biased systems. The grain size dependence of the shows the expected peak due to a competition between the superparamagnetic nature of small grains and reduction in the statistical imbalance in the number of interfacial spins for larger grain sizes. Our simulations confirm the existence of single antiferromagnetic domains within each grain. The model gives insights into the physical origin of exchange bias and provides a route to developing optimised nanoscale antiferromagnetic spintronic devices.