Abstract: Based on the spectral analysis theory and the turbulent physical spectral tensor model, a stochastic harmonic function method is proposed for the simulation of 3-directional (3D) fluctuating wind fields. The physical spectral tensor model proposed in the IEC code is firstly introduced based on the physical properties and the rapid distortion theory (RDT) of the 3D atmospheric boundary layer turbulent field. Then, by introducing the stochastic harmonic function representation, the 3D fluctuating field is expressed as the superposition of a series of harmonic functions modulated by the stochastic wavenumbers and phase angles. A two-step acceptance-rejection scheme is proposed to determine the representative point set of the stochastic wavenumbers. The zero-phase evolution time model is adopted to express the high-dimensional stochastic phase angles as the stochastic functions of 4-dimentional random variables. For the purpose of verification, the numerical examples of simulating the 3D fluctuating wind field for a bridge tower are performed. It is demonstrated that the proposed method can simulate the 3D fluctuating wind field corresponding to different turbulence modes with a small number of harmonic terms and fair accuracy.