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.
ndsolve
(not to be confused withNDSolve
) and\[Sigma]gr
. Please clarify. $\endgroup$