The blade-tower interaction of upwind horizontal axis wind turbines has become important to aerodynamic loading as the systems become larger. However, there are not enough studies describing these phenomena. To investigate this interaction, we performed numerical simulations for uniform, yawed, wind shear flow conditions, and various tower cases using the nonlinear vortex correction method with time-marching free wake. At 5 m/s, the change in the normal force coefficient is approximately 10% of the average. The blade root region has a larger azimuth range of the interaction and a bigger change in aerodynamic loading. The blade-tower interaction decreases as the yaw error and wind shear exponent increases. The interaction due to tower radius variations is higher than that due to tower clearance variations. With regard to stochastic load, the blade-tower interaction may affect the total fatigue load at low wind speed and in a more unstable atmospheric condition.