Due to neutron irradiation, solid breeder blankets are subjected to complex thermo-mechanical conditions. Within one breeder unit, the ceramic breeder bed is composed of spherical-shaped lithium orthosilicate pebbles, and as a type of granular material, it exhibits strong coupling between temperature and stress fields. In this paper, we study these thermo-mechanical problems by developing a thermal discrete element method (Thermal-DEM). This proposed simulation tool models each individual ceramic pebble as one element and considers grain-scale thermo-mechanical interactions between elements. A small section of solid breeder pebble bed in a helium-cooled pebble bed (HCPB) is modelled using thousands of individual pebbles and subjected to volumetric heating profiles calculated from neutronics under ITER-relevant conditions. We consider heat transfer at the grain scale between pebbles through both solid-to-solid contacts and the interstitial gas phase, and we calculate stresses arising from thermal expansion of pebbles. The overall effective conductivity of the bed depends on the resulting compressive stress state during the neutronic heating. The Thermal-DEM method proposed in this study provides access to the grain-scale information, which is beneficial for HCPB design and breeder material optimization, and a better understanding of overall thermo-mechanical responses of the breeder units under fusion-relevant conditions.