Due to the agile maneuverability, unmanned aerial vehicles (UAVs) have shown great promise for on-demand communications in the next-generation wireless networks. Considering the massive multiple-input multiple-output (MIMO) configuration, this paper proposes a novel three-dimensional (3D) beam domain channel model (BDCM) for UAV communications. Through dividing the large antenna array into several sub-arrays and classifying multipath components as near-field and far-field components, the proposed BDCM takes the spherical wave front (SWF) and array non-stationarity into account. Channel statistical properties including spatial-temporal-frequency correlation function (STF-CF), root-mean-squared (RMS) Doppler spread, beam spread, channel matrix collinearity (CMC), and stationary time interval are derived and simulated for the proposed BDCM. Influences of SFW and non-stationary properties on the statistical properties and system performance are analyzed. Simulation results show that, compared with the equivalent geometry-based stochastic model (GBSM), the proposed BDCM has better temporal correlation, while BDCM and GBSM are equivalent in the system performance evaluation. Furthermore, the performance of the proposed BDCM is evaluated in terms of accuracy, complexity, and pervasiveness. The results show that the proposed BDCM can represent massive MIMO channel properties accurately with low complexity and good compatibility.