piecewise - Solving a discontinuous differential-algebraic equation system for plasticity behaviour

I need to solve a discontinuous equation which is typical in theory of plasticity. For a simple case I get the following equation system (reformulated for numerical implementation):

$$\begin{align*} s(t) &= \frac{\sigma(t)}{C_1} + s_{ep}\\ s_{ep}'(t) &= \begin{cases}\frac{C_1}{C_1+C_2}s'(t) & \text{for } |\sigma(t)-C_2 s_{ep}(t)| \ge \sigma_{gr} \land \sigma(t)s'(t)>0\\0 & \text{otherwise}\end{cases} \end{align*}$$

with "zero" initial conditions. I'd like to get the solution for $s(t)$ for given parameters $C_1$, $C_2$, $\sigma_{gr}$ and a known function $\sigma(t)$. I assumed: $\sigma(t) = 40000\sin(0.02t)$, $C_1=80000$, $C_2 = 20000$, $\sigma_{gr} = 15000$. This should give a hysteresis loop on a plane $\sigma(t)-s(t)$.

So in Mathematica I tried to use automatic discontinuity handling by defining the second equation using a Piecewise function:

σ[t_] := 40000*Sin[0.02*t];
eq1 = s[t] == σ[t]/C1 + sep[t];
eq2 = sep'[t] == Piecewise[{{C1/(C1 + C2)*s'[t], (σ[t]*s'[t] > 0) && ((σ[t] - C2*sep[t] >=  σgr) || (σ[t] - C2*sep[t]<=-σgr))}}, 0];
eqSys := {eq1, eq2, s[0] == 0, sep[0] == 0};
ndsolve=NDSolve[eqSys, {s[t], sep[t]}, {t, 0, 1000}]
disp[t_] := Evaluate[s[t] /. ndsolve];
sTab = Table[disp[t][[1]], {t, 0, 1000, 1}];
σTab = Table[σ[t], {t, 0, 1000, 1}];
ListPlot[Transpose[{sTab, σTab}], PlotRange -> All, GridLines -> Automatic]

Unfortunately I get:

NDSolve::tddisc: NDSolve cannot do a discontinuity replacement for event surfaces that depend only on time. >>

and the results are incomplete or the algorithm crashes. I also tried using WhenEvent with "DiscontinuitySignature" but with no success. This approach gives good results only for a linear monotonic function of $\sigma$, e.g. $\sigma(t) = 50t$.

I wrote a module to solve this using a simple first order Runge-Kutta so I obtained the solution but this is only a simple model. I'm sure Mathematica can solve this with its build-in methods. That would really save me a lot of work writing my own procedures.

Answer

Edit:
Using a helper function fh will result in no messages and no need to set extra options.

σ[t_] := 40000 Sin[0.02 t]

C1 = 80000;
C2 = 20000;
σgr = 15000;

fh[t_?NumericQ, x_, y_] := Piecewise[{{C1/(C1 + C2)*y, 
  (σ[t]*y > 0) && ((σ[t] - C2*x >= σgr) || (σ[t] - C2*x <= -σgr))}}, 0]

sol = NDSolve[{s[t] == σ[t]/C1 + sep[t], sep'[t] == fh[t, sep[t], s'[t]], 
        s[0] == 0, sep[0] == 0}, {s[t], sep[t]}, {t, 0, 1000}];


s[t_] = s[t] /. sol // First;
ParametricPlot[{s[t], σ[t]}, {t, 0, 10^3}, PlotRange -> All, 
  AspectRatio -> Full, GridLines -> Automatic]