Solving a discontinuous differential-algebraic equation system for plasticity behaviour

Question

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):

$s(t) = \frac{\sigma(t)}{C_1} + s_{ep}$

$s_{ep}'(t) = \left\{\begin{array}{lll}\frac{C_1}{C_1+C_2}s'(t) & \textrm{for}& |\sigma(t)-C_2 s_{ep}(t)| \ge \sigma_{gr} \land \sigma(t)s'(t)>0\\0 & & otherwise\end{array}\right.$

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.02*t)$, $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) = 50*t$.

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 1

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]

Answer 2

I have been unsuccessful at getting the automatic discontinuity handling to work. (I get the errors "Failure to project onto the discontinuity surface when computing Filippov continuation".) But manually handling it with WhenEvent works, although it complains about slow convergence to the even locations. Perhaps the discontinuity conditions are too complicated for smooth handling. I don't know.

σ[t_] := 40000 Sin[1/50 t];
Block[{s, C1 = 80000, C2 = 20000, σgr = 15000},
 eq1 = s[t] == σ[t]/C1 + sep[t];
 s[t_] := σ[t]/C1 + sep[t];
 eq2 = sep'[t] == sepprime[t] C1/(C1 + C2)*s'[t];
 events = Simplify[
    Solve[
     sep'[t] == Piecewise[{{C1/(C1 + C2)*s'[t], (σ[t]*s'[t] > 0) &&
         ((σ[t] - C2*sep[t] >= σgr) || (σ[t] - C2*sep[t] <= -σgr))}}, 0],
     sep'[t]],
    (sep[t] | σ'[t]) ∈ Reals] /. 
   HoldPattern[{sep'[t] -> ConditionalExpression[val_, cond_]}] :> 
    WhenEvent[cond, sepprime[t] -> If[val === 0, 0, 1]];
 ]
eqSys = {eq1, eq2, events, s[0] == 0, sep[0] == 0, sepprime[0] == 0};

sol = Quiet[
  NDSolve[eqSys, {s, sep}, {t, 0, 1000}, DiscreteVariables -> {sepprime}],
  NDSolve::evcvmit];

OP's plot:

ParametricPlot @@ {{s[t], σ[t]} /. First[sol], 
  Flatten[{t, sep["Domain"] /. First[sol]}], AspectRatio -> 2/3}

The DE can also be integrated by turning off the discontinuity handling, which is effectively similar to Karsten 7's solution.

sol = NDSolve[eqSys, {s, sep}, {t, 0, 1000}, Method -> {"DiscontinuityProcessing" -> False}]

This is effectively what the OP's code does, after the NDSolve::tddisc warning message, but it does not crash in V10.0.1. In V9, it one needs to set MaxSteps -> 60000 or StartingStepSize -> 0.1 for the integration to complete.