In this work, a multiscale data-driven (MSDD) methodology is proposed to efficiently simulate the thermomechanical behavior of densely packed, confined granular materials subjected to thermal expansion effects [1,2]. The multiscale strategy employed follows the hierarchical approach, where the macroscale is handled using a continuous computational model, particularly the Finite Volume Method (FVM), while the microscale response is obtained from Representative Volume Elements (RVEs) with the Discrete Element Method (DEM) [3]. To significantly reduce the computational cost of the analyses, the microscale DEM computations are not performed online, i.e. simultaneously with the macroscale solver, as commonly done in the standard hierarchical multiscale methodology. Instead, they are performed in advance to create a comprehensive database of RVE solutions under different initial conditions and thermal strains. This dataset is then used to train an Artificial Neural Network (ANN), which serves as a surrogate model for the macroscale solver. The MSDD approach is validated against pure DEM solutions of two-dimensional enclosed granular domains with distinct thermal boundary conditions. Remarkably, we demonstrate that with only three input parameters, namely porosity, fabric index, and thermal strain, the surrogate model can predict the evolution of the microstructure, as well as the updated tensors of effective thermal conductivity and Cauchy stress for the granular assembly. This allows for a generally accurate simulation of transient thermomechanical analyses at a drastically lower computational cost than the pure DEM approach.