An explicit analytic solution for thermal stresses in an infinite thermoelastic medium with two circular cylindrical holes of different sizes kept at different constant temperatures, under steady-state heat flux is presented. The solution is obtained by using the most general representation of a biharmonic function in bipolar coordinates. The stress field is decomposed into the sum of a particular stress field induced by the steady-state temperature distribution and an auxiliary isothermal stress field required to satisfy the boundary conditions on the holes. The variations of the stress concentration factor on the surface of the holes are determined for varying geometry of the holes. The concept of the conservation integrals Jk, M and L is extended to steady state thermoelasticity and the integrals are proved to be path-independent. These integrals are calculated on closed contours encircling one or both holes. The geometries of a hole in a half-space and an eccentric annular cylinder are considered as particular cases.