The present work focus on the interlaminar stress formulation based on equilibrium equations, applied in a new triangular multi-layer nonlinear shell finite element suitable for simulation with large displacements and rotations. The shell element used is an adaptation to a multilayer situation of the T6-3iKL element developed in [1], a kinematical model with properties from Kirchhoff-Love shell theory, considering the shell director across the layers as constant and the Rotation-continuity between adjacent elements and allowing multiple branches connections in the mesh. Additionally, the element considers strains along the edge of the element, as an extrapolation of the work developed in [2] and [3], and an additional bubble node as mean to guarantee the necessary degree of the interpolation to analyse the stress distribution along the thickness. The element is developed allowing for implementation of different material constitutive equations (saint-vennant, Neo-hookean and anisotropic materials). The model developed in this article is numerically implemented and results are compared to different references in multiple examples, showing the consistency and robustness of the formulation. Regarding the interlaminar stress, a crank-nicolson equilibrium based scheme is developed, this being an extremelly simple approach, capable of representing complex stress distributions and stablishing the foundations for future work on delamination analyses. It is believed that the multilayer extension with the desirable properties of no necessity of artificial penalty calibration, simple kinematic, a relatively small number of DOFs, geometric exact, possibility to use 3D material constitutive models, easily connected with multiple branched shells and beams, and including possibly the most simple consideration of multilayers and interlaminar stress calculation, create a simple yet powerfull shell element for interlaminar stress analyses.