Solving a Nonlinear Complementary Problem (plasticity)

I think the behavior of WhenEvent encountered by OP is a bug. Anyway, here's a working WhenEvent based solution:

sigma[t_] = Sin[t];
sigmay = 0.5;
E0 = 1;
g[t_, epsi_] = sigmay - Abs[sigma[t] - E0*epsi];
phi[t_, epsi_, dotepsi_] = (sigma[t] - E0 epsi) dotepsi;

events = {WhenEvent[g[t, epsi[t]] < phi[t, epsi[t], epsi'[t]], coef[t] -> 1], 
   WhenEvent[phi[t, epsi[t], epsi'[t]] < g[t, epsi[t]], coef[t] -> 0]};

epsisol = First@
  NDSolveValue[{g[t, epsi[t]] coef[t] + phi[t, epsi[t], epsi'[t]] (1 - coef[t]) == 0, 
    epsi[0] == 0, coef[0] == 0, events}, {epsi, coef}, {t, 0, 100}, 
   DiscreteVariables -> coef, SolveDelayed -> True]

Plot[epsisol[t], {t, 0, 15}]

A trick to obtain the full result.

sigma[t_] := Sin[t];
sigmay = 0.5;
E0 = 1;
tmax = Pi;
g[t_?NumericQ, epsi_] := sigmay - Abs[sigma[t] - E0*epsi]
phi[t_?NumericQ, epsi_, dotepsi_] := (sigma[t] - E0*epsi)*dotepsi
tmax = Pi;
tmin = 0;
epsisolant = sigma[tmin];
GR = {};

While[tmax < 100, 
  epsisol = NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, epsi[tmin] == epsisolant}, epsi, {t, tmin, tmax}, Method -> {"EquationSimplification" -> "Residual"}];
  AppendTo[GR, Plot[epsisol[t], {t, tmin, tmax}]];
  epsisolant = epsisol[tmax];
  tmin = tmax;
  tmax += Pi/4
]

Show[GR, PlotRange -> All]

You input is to my knowledge applied to it correctly. Well done.

But this is a discretized attempt to solve the problem.

sigma[t_] := Sin[t];
sigmay = 0.5;
E0 = 1;
g[t_?NumericQ, epsi_] := sigmay - Abs[sigma[t] - E0*epsi]
phi[t_?NumericQ, epsi_, dotepsi_] := (sigma[t] - E0*epsi)*dotepsi
epsisol = 
 NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, 
   epsi[0] == 0}, epsi, {t, 10^-13, 100}]

The second message opens a page ndsolve::ndcf with the direct invitation to contact the technical support of Wolfram Inc.

I found that the domain depends with a rapid jump on the starting time at a little bit more than 10^-13 for example a quarter I reproduce Your results and around that less again. Might be this is a match for the domain length 4.71. That can even be got again at higher starting times as 0.0001 or so.

My output is:

Plot[epsisol[t], {t, 0.005, 4.71}, PlotRange -> Full]

From that on I agree with [@cesareo]5 it might go on delayed quasi-periodic. This might already be chaotic not only in the starting time but in the period. The rise and fall may be characteristic. Somehow this is similar to a sawtooth. Therefore and because the switch function suggests it, I make the solution ideation that this might be solved with Fourier or Laplace methodologies for more domain. This will only work in approximation.

But the curious idea changed my plans: make the domain smaller arbitrarly:

epsisol = 
 NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, 
   epsi[0] == 0}, epsi, {t, 10^-13, 10}]

Plot[epsisol[t], {t, 0.005, 10}, PlotRange -> Full]

Hope that helps. This is done with V12.0.0 on iMac Catalina.

This can be solved up to 10.99639 if the Method -> {"EquationSimplification" -> "Residual"} is used. The message remains: ndcf. The repeated convergence test does not accept the rapid stagnation of the growth of the solution at -0.5. But it suffices for the full period of the graph. Perhaps the treatment as a differential-algebraic equation.

Seems that a better match for sigmay and sigma give a longer domain in the capabilities for of-the-shelf differential-algebraic methods. Maybe this is on the other-hand a question designed for th fail of the adaptibility of the repeated convergence test.

I made a comparison between both solution, mine and from @xzczd.

Plot[{epsisol[t], epsisolu[t], 
  1.25 10^7 (epsisol[t] - epsisolu[t])}, {t, 0.00001, 10.99}, 
 PlotRange -> Full]

Despite both solutions look at first sight very similiar they are different.

Plot[{epsisol[t] - .5, epsisolu[t] - .5}, {t, 1.8, 3.14}, 
 PlotRange -> Full, PlotLegends -> "Expressions"]

Mine stays a little bit, an order of one magnitude further away from the limiting 0.5. This is even bigger for the negative border and bigger at the second constant interval. Then my solution goes over to fail. The even very small error oscillates up and the end is the test fails.

With InterpolationOrder->All the oscillation get much smaller and more repetitive:

But the domain is not larger.

To each Accuracy 9,10,11,... there is interval near zero for which the integration is successful.

epsisol = NDSolveValue[{Min[gi[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, epsi[0] == 0}, epsi, {t, 10^-10.1295, 11}, Method -> {"EquationSimplification" -> "Residual"}, InterpolationOrder -> All, AccuracyGoal -> 10]

Plot[{epsisol[t], epsisolu[t], 
  1.25 10^7 (epsisol[t] - epsisolu[t])}, {t, 0.00001, 7.85}, 
 PlotRange -> Full]

For Accuracy 11 the domain has a much large interval for which my solution gets much closer to the reference solution and the oscillation get tamed. At -0.5 mine is better than that of the competitor. But the oscillations remain still order 10^-7.

Fast and dirty as Mathematica built-ins are these days. The behavior is a clear hint that Mathematica uses StiffnessSwitching internally for the calculation of the solution.

ListLinePlot@
 Quiet@Table[(epsisol = 
      NDSolveValue[{Min[gi[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 
         0, epsi[0] == 0}, epsi, {t, 10^expon, 11}, 
       Method -> {"EquationSimplification" -> "Residual"}, 
       InterpolationOrder -> All, AccuracyGoal -> 13])[[1, 1, 
     2]], {expon, -5, -16, -.01}]

There are many possible start values for Accuracy 12. The result is still stiffness switching wildly but the exactness increases strongly.