Abstract: The upper atmosphere density is the primary source of prediction error for the aerodynamic drag of spacecraft in ultra-low orbit. The atmospheric density is usually provided by semi-empirical atmospheric models, and its uncertainty directly affects the accuracy of the drag. Taking the GOCE satellite as the research object, the aerodynamic drag analysis was carried out using the test particle Monte Carlo(TPMC) method. A calculation scheme for the aerodynamic drag of the GOCE satellite in orbit was proposed. By comparing with the in-orbit measured data of drag, the accuracy of different atmospheric models was quantitatively evaluated. The influence of spatial variables such as orbit altitude, latitude, and longitude, as well as time variables such as seasonal alternation and day-night cycle, on the aerodynamic drag was systematically studied. The results show that compared with the in-orbit measured data of drag of the GOCE satellite, the maximum relative deviation of the drag predicted by the five atmospheric models does not exceed 40%, and the average relative deviation does not exceed 20%. Among them, the average deviation of USSA-1976 is the largest, and that of NRLMSIS-00 is the smallest. In the prediction of aerodynamic drag of ultra-low orbit satellites, it performs the most accurately. The aerodynamic drag of ultra-low orbit satellites decreases exponentially with the increase in altitude. With the enhancement of solar activity, the aerodynamic drag of the satellite gradually increases. The higher the orbit altitude, the more obvious the influence of solar activity. The influence of seasonal alternation and day-night cycle on the aerodynamic drag of ultra-low orbit satellites is reflected through the solar radiation intensity. The greater the radiation intensity at the satellite’s latitude and longitude position and the season, the greater the drag. When the spacecraft moves east-west, the drag fluctuates more than when it moves north-south. The relative change of drag within a day exceeds 30%.