I would like to solve the following: given $t\mapsto\sigma(t)$ and $E>0$, $\sigma_y>0$, find $\epsilon$ such that:
$$\left\lbrace\begin{array}{l}g(t,\epsilon)\geq 0,\\ \phi(t,\epsilon,\epsilon')\geq0, \\g(t,\epsilon)\phi(t,\epsilon,\epsilon')=0,\\ \epsilon(0)=0. \end{array}\right.$$
with $g(t,\epsilon(t))=\sigma_y - |\sigma(t) - E\epsilon(t)|$ and $\phi(t,\epsilon(t),\epsilon'(t)) = (\sigma(t) - E\epsilon(t))\epsilon'(t)$.
Example:
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, 0, 100}] (* NDSolveValue::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. NDSolveValue::ndcf: Repeated convergence test failure at t == 1.5717016618338167`; unable to continue. *) Attempt with WhenEvent:
events = {WhenEvent[g[t, epsi[t]] == 0, coef[t] -> 1], WhenEvent[phi[t, epsi[t], epsi'[t]] == 0, 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, 10}, DiscreteVariables -> {coef}] (* NDSolveValue::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. ... then integration stops at t = 1.54 with no further error *) Context and details
I am interested in plasticity, in particular the computation of the displacement $\varepsilon$ induced by an applied force $\sigma$ to a system composed of a slider and a spring in parallel, such as:
This problem is addressed in Solving a discontinuous differential-algebraic equation system for plasticity behaviour ($C_2$ is $H$ and $C_1$ is $\infty$) --- I believe there is a slight mistake in the equations but it still works after correction.
However, with both answers, I did not manage to adapt the code so that it would work with two such devices in series. Also, I wanted to derive the equation from "standard" plasticity theory:
- a plasticity criterion $f =|\sigma -X| - \sigma_y \leq 0$ ($g=-f$ above, and $X=E \epsilon$ is the stress in the spring);
- the positiveness of dissipation: $\phi = \epsilon' (\sigma_y -E\epsilon)\geq 0$
- an orthogonality condition implying that energy is dissipated iff there is plasticity ($f=0$): $f\times \phi = 0$.
This is often written altogether: $$ 0\leq (-f) \perp \phi \geq 0$$
Such kinds of formulation are also found in intermittent contact dynamics: the reaction force is always non-negative, and can be non-zero only if there is contact, i.e. when the distance is zero. Conversely, if the distance is non-negative, the reaction for can only be zero.
Such problems are numerically challenging, even though there are dedicated numerical methods. Even the formulation involving the derivative $\epsilon'$ is wobbly, because $\epsilon$ is not differentiable everywhere (just like the velocity of a bouncing ball is not defined at impact times).
Anyway, WhenEvent works very fine for bouncing balls with few contacts, so I would have expected WhenEvent to be efficient here.
