# 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, sep == 0};
ndsolve=NDSolve[eqSys, {s[t], sep[t]}, {t, 0, 1000}]
disp[t_] := Evaluate[s[t] /. ndsolve];
sTab = Table[disp[t][], {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.

• Your code above has two undefined quantities, ndsolve (not to be confused with NDSolve) and \[Sigma]gr. Please clarify. – bbgodfrey Dec 23 '14 at 14:41
• @Karsten7. It does! I get the same s(t) when I use my "Euler method" module. It should be hysteresis loop when you plot a function s(sigma). – K.J. Dec 23 '14 at 14:45
• For the record, I believe there is a mistake in the equations: it should be $(\sigma(t)-C_2 s_{ep}(t))s'(t)>0$ instead of $\sigma(t) s'(t)>0$. – anderstood Oct 21 at 15:52
• @anderstood You are right. I'm not going to edit the question though as the mistake is not related to solving this kind of problem in Mathematica. – K.J. Nov 20 at 15:00

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, sep == 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] • What changes did you make? When I copy-paste your code to my Mathematica (ver. 9.0) I get such a graph: link and an error: NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 19.964613731417995. This is what I obtain using my module:link. – K.J. Dec 23 '14 at 15:18
• @K.J., I can not reproduce the crash in 10.0.2 seems fixed. – user21 Dec 23 '14 at 16:28
• @Karsten 7. Yes, with increased MaxSteps I get correct results though the first error is still present as you've written. It appears that is the case. However the calculation takes a while. I hope it won't crash for a more complicated system of diff. equations. I'll check that. – K.J. Dec 23 '14 at 16:28
• @K.J. Please see my edit. With this approach no messages are raised and no extra options need to be set. You can find an explanation for this strategy here. – Karsten 7. Dec 23 '14 at 17:41
• @Karsten7. Thanks a lot. That works perfectly indeed. – K.J. Dec 24 '14 at 10:30

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 event 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, sep == 0, sepprime == 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, one needs to set MaxSteps -> 60000 or StartingStepSize -> 0.1` for the integration to complete.