A relatively simple two-dimensional multilayered shell model is presented for predicting both global quantities and stress distributions across the thickness of multilayered thick shells, that is based on a third-order zig-zag approach. As for any zig-zag model, the layerwise kinematics is accounted for, with the stress continuity conditions at interfaces met a priori. Moreover, the shell model satisfies the zero transverse shear stress conditions at the upper and lower free surfaces of the shell, irrespective of the lay-up. By changing the parameters in the displacement model, some higher order shell models are obtained as particular cases. Although it potentially has a wide range of validity, application is limited to cylindrical shell panels in cylindrical bending, a lot of solutions of two-dimensional models based on rather different simplyfying assumptions and the exact three-dimensional elasticity solution being available for comparisons for this benchmark problem. The numerical investigation performed by the present shell model and by the shell models derived from it illustrates the effects of transverse shear modeling and the range of applicability of the simplyfying assumptions introduced. The implications of retaining only selected terms depending on the radius-to-thickness ratio are focused by comparing the present solutions to the exact one and to other two-dimensional solutions in literature based on rather different simplyfying assumptions.