Multilayered plate structures, such as composite laminates, are regarded as high-performance material solutions in many engineering fields, e.g. aerospace, naval, and automotive industries. Providing high-fidelity simulations at low computational cost is regarded as a rewarding objective to accelerate the design of new allowables. To address this challenge, the present work extends the multiscale 2D+ methodology introduced in [1] to the nonlinear analysis of multilayered flat shells under bending within an industrially relevant finite element environment. The proposed framework exploits dimensional reduction at both the macro-and meso-scales, significantly decreasing computational cost while retaining the heterogeneous through-thickness response. It relies on a computational homogenization strategy in which flat shell structures with degenerated kinematics are employed at the macro-scale, coupled with a one-dimensional Representative Volume Element (RVE) across the thickness at the meso-scale. Representative case studies employing a 2D+ flat shell formulation implemented in Abaqus demonstrate: (i) stress distribution predictions comparable in accuracy to those obtained from full three-dimensional simulations, (ii) computational efficiency equivalent to that of conventional two-dimensional approaches based on Equivalent Single Layer (ESL) theories, (iii) a fully non-intrusive implementation restricted to the element level via user-defined elements (UEL), and (iv) the ability to naturally capture material nonlinear behavior.