I am modelling simple periodic shear on a C3D8 element, finite strain, Neo Hooke material with Prony series. I have noticed an evident discrepancy in stress prediction (factor 10^5) when g_1 (one branch) tends to 0.999, between Abaqus and analytical solution (which is matched by my python prediction following Simo's algorithms). Please, find attach my .inp file. 

Since it is documented that Abaqus should follow Simo's algorithms I would expect a good match with the analytical solution (as Simo's algorithms in my python script do) even for g_1=0.999. I cannot find an explanation to this problem in literature. Any suggestionAbaqus