Up to date, accurate prediction of interlaminar stresses is still a challenging issue for two-node beam elements. The postprocessing approaches by integrating the three-dimensional equilibrium equation have to be used to obtain improved transverse shear stresses, whereas the equilibrium approach requires the first-order derivatives of in-plane stresses. In-plane stresses within two-node beam element are constant, so the first-derivatives of in-plane stresses are close to zero. Thus, two-node beam elements encounter difficulties for accurate prediction of transverse shear stresses by the constitutive equation or the equilibrium equation, so a robust two-node beam element is expected. A two-node beam element in terms of the global higher-order zig-zag model is firstly developed by employing the three-field Hu-Washizu mixed variational principle. By studying the effects of different boundary conditions, stacking sequence and loading on interlaminar stresses of multilayered composite beams, it is shown that the proposed two-node beam element yields more accurate results with lesser computational cost compared to various higher-order models. It is more important that accurate transverse shear stress has active impact on displacements and in-plane stresses of multilayered composite beams.