it seems that I’ve encountered a bug. The problem is that the results from a multistep analysis (2 identical static steps without NLGEOM) involving shells and OP=NEW used for a boundary condition applied to translational and rotational DOFs (the latter seem to be the source of the problem since this issue doesn’t occur with solid elements) are incorrect in the second step. Abaqus gives correct solution for the same input file.
Good idea but if I add this boundary condition before the first step (keeping the remaining definitions within the steps because OP=NEW means that the BC has to be redefined) the results are still the same. So there seems to be a different reason for this bug. I think that I will report it in the GitHub repository.