Channel estimation has long been deemed as one of the most critical problems in three-dimensional (3D) massive multiple-input multiple-output (MIMO), which is recognized as the leading technology that enables 3D spatial signal processing in the fifth-generation (5G) wireless communications and beyond. Recently, by exploring the angular channel model and tensor decompositions, the accuracy of single-user channel estimation for 3D massive MIMO communications has been significantly improved given a limited number of pilot signals. However, these existing approaches cannot be straightforwardly extended to the multi-user channel estimation task, where the base station (BS) aims at acquiring the channels of multiple users at the same time. The difficulty is that the coupling among multiple users' channels makes the channel estimation deviate from widely-used tensor decompositions. It gives a non-standard tensor decomposition format that has not been well tackled. To overcome this challenge, besides directly fitting the new tensor model for channel estimation to the wireless data via block coordinate descent (BCD) method, which is prone to the overfitting of noises or requires regularization parameter tuning, we further propose a novel tuning-free channel estimation algorithm that can automatically control the channel model complexity and thus effectively avoid the overfitting. Numerical results are presented to demonstrate the excellent performance of the proposed algorithm in terms of both estimation accuracy and overfitting avoidance.