A bolted joint is frequently used under thermal load in practical applications, such as pressure vessels, internal combustion engines, etc. In order to accurately evaluate the thermal stresses thus produced, the effects of thermal contact resistance at the interface and the heat flow through small gaps, which exist around the objective bolted joint, must be taken into account. In this paper, a numerical approach with high computation efficiency is proposed, where empirical equations for thermal contact coefficient and apparent thermal contact coefficient are incorporated into commercial engineering software. By conducting systematic three-dimensional FE analyses, it is quantitatively elucidated how the supplied heat flows through each part of a bolted joint and how the bolt stress varies with time. Bolted joints made of the materials with low thermal conductivity exhibit specific behaviors on the heat flow pattern around the bolted joint and the variations of axial bolt stress and bolt bending stress.