Refractory ceramics, widely used in high-temperature industrial applications, exhibit highly heterogeneous microstructures composed of multi-phase compositions, aggregates of various sizes, and different bonding systems. This complexity results from their exposure to extreme environments, where they must withstand thermal shocks, thermal cycling, severe thermal gradients, mechanical loads, corrosion, and abrasion. A critical property for ensuring durability under these conditions is thermal shock resistance, which depends on key macroscopic thermomechanical properties such as stiffness, fracture toughness, Poisson's ratio, thermal expansion coefficient, and thermal conductivity. These properties are significantly influenced by the presence of microcracks, which may be intentionally introduced to enhance energy dissipation during fracture. Among refractory materials, aluminum titanate is particularly relevant due to its unique crystal structure conducting to strongly anisotropic thermal expansion behavior at the grain level. During high-temperature operations, mismatches in the coefficient of thermal expansion (CTE) between grains induce spontaneous microcracking, leading to a quasi-brittle, non-linear mechanical response. This non-linearity is often associated with material toughening, enhancing fracture resistance and improving thermal shock performance. Traditional approaches for modeling microcrack behavior rely on continuum-based methods such as the finite element method with cohesive zone models, continuum damage mechanics, phase-field models, or extended finite element methods. However, these techniques face fundamental challenges in handling multi-fracture problems, such as crack initiation, propagation, branching, multi-cracking, and crack closure. An alternative approach is the discrete element method, which inherently manages discontinuities and has been successfully applied to rocks and concrete-like materials. Despite its advantages, few discrete element method studies incorporate thermal effects at the microstructural scale, which are key to predicting microcrack networks and their impact on macroscopic properties. To address this limitation, this study employs a novel discrete element method-based approach, using a bonded particle model and integrating anisotropic thermal expansion, thermomechanical coupling, and crack closure mechanisms under periodic boundary conditions to enable a multiscale analysis. This methodology extends