the distributed normal force (n) in a shell cross section is the force per unit of length. In order to convert it to a stress (sigma) one needs to divide it by the shell thickness t: sigma = n / t.
Distributed bending moments (m) give rise to a linear stress distribution over the shell thickness such that sigma = m * 12 / (t^3)*z. Where z is the distance from the cross section center. The maximum stress at the upper and lower boundary of the shell (z = ±t/2) thus comes out as sigma = m * 6 / t^2.
The principal directions of the moments and normal forces do not necessarily coincide. Therefore one can not simply add the contributions from e.g. n1 and m1. Instead one needs to start from the distributed forces and moments in the local element coordinate system, calculate the stress state at the position in question and then calculate the principal stresses.
An alternative would be to use Karamba3D 2.0.0 WIP (see here), specify a material with a Rankine yield criterion, set the tensile and compressive strength to ‘1’ and calculate the utilization of the shell elements.