The temperature distribution throughout the LAGEOS satellites is simulated numerically with the objective to determine the resulting thermal force. The different elements and materials comprising the spacecraft, with their energy transfer, have been modeled with unprecedented detail. The radiation inputs on the satellites are direct solar (eclipse modulated), Earth albedo, and Earth infrared radiations. For each satellite the lifetime temperature (behavior) of 2133 nodes is computed. On the basis of this distribution, individual forces and the net instantaneous accelerations are obtained. Simulations yield typical temperature variations ranging between 30 and 100 K for different elements and materials, whereas the net instantaneous accelerations are on the order of 70 pm s-2, in good agreement with previous results. Simulations also show the importance of the consideration of a proper orientation of the satellite: LOSSAM yields acceleration differences of up to three times the acceleration obtained with a constant spin axis orientation. The temperature of the four germanium retroreflectors deviates up to 70 and 100 K with respect to their silica counterparts for LAGEOS I and II, respectively. This generates thermal acceleration differences of several pm s-2, up to 25% of the postulated difference in reflectivity between hemispheres. Two factors play a major role: the spin rate and the Sun aspect angle with respect to the spin axis. On the basis of the latter, two characteristic periods can be distinguished: a rapid spin, slow drift period (until 13 and 8 years after launch for LAGEOS I and II, respectively) and a slow spin, rapid wobbling afterward. The acceleration results will be used in a refined orbit computation in a subsequent investigation.