Lösung eines nichtlinearen komplementären Problems (Plastizität)

Question

Ich möchte folgendes lösen: gegeben $t\mapsto\sigma(t)$ und $E>0$, $\sigma_y>0$, finden $\epsilon$ so dass:

$$\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.$$

mit $g(t,\epsilon(t))=\sigma_y - |\sigma(t) - E\epsilon(t)|$ und $\phi(t,\epsilon(t),\epsilon'(t)) = (\sigma(t) - E\epsilon(t))\epsilon'(t)$.

Beispiel:

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. *)

Versuch mit 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 *)

---

Kontext und Details

Ich interessiere mich für Plastizität, insbesondere die Berechnung der Verschiebung $\varepsilon$ induziert durch eine aufgebrachte Kraft $\sigma$ zu einem System, das aus einem Schieber und einer Feder parallel besteht, wie z.

Dieses Problem wird in der Lösung eines diskontinuierlichen differentiell-algebraischen Gleichungssystems für das Plastizitätsverhalten ($C_2$ ist $H$ und $C_1$ ist $\infty$) --- Ich glaube, es gibt einen kleinen Fehler in den Gleichungen, aber es funktioniert immer noch nach der Korrektur.

Bei beiden Antworten gelang es mir jedoch nicht, den Code so anzupassen, dass er mit zwei solchen Geräten in Reihe funktioniert. Außerdem wollte ich die Gleichung aus der "Standard" -Plastizitätstheorie ableiten:

• ein Plastizitätskriterium $f =|\sigma -X| - \sigma_y \leq 0$ (($g=-f$ oben und $X=E \epsilon$ ist der Stress im Frühjahr);
• die positive Wirkung der Dissipation: $\phi = \epsilon' (\sigma_y -E\epsilon)\geq 0$
• eine Orthogonalitätsbedingung, die impliziert, dass Energie bei Plastizität abgeführt wird ($f=0$): $f\times \phi = 0$.

Dies wird oft zusammen geschrieben: $$ 0\leq (-f) \perp \phi \geq 0$$

Solche Formulierungen finden sich auch in der intermittierenden Kontaktdynamik: Die Reaktionskraft ist immer nicht negativ und kann nur bei Kontakt ungleich Null sein, dh wenn der Abstand Null ist. Wenn umgekehrt der Abstand nicht negativ ist, kann die Reaktion für nur Null sein.

Solche Probleme sind numerisch herausfordernd, obwohl es spezielle numerische Methoden gibt. Sogar die Formulierung mit dem Derivat$\epsilon'$ ist wackelig, weil $\epsilon$ ist nicht überall differenzierbar (genau wie die Geschwindigkeit eines springenden Balls zu Aufprallzeiten nicht definiert ist).

Wie auch immer, WhenEventfunktioniert sehr gut zum Hüpfen von Bällen mit wenigen Kontakten, daher hätte ich erwartet, dass ich WhenEventhier effizient bin .

Answer 1

Ich denke, das Verhalten WhenEventvon OP ist ein Fehler. Hier ist eine funktionierende WhenEventLösung:

sigma[t_] = Sin[t];
sigmay = 0.5;
E0 = 1;
g[t_, epsi_] = sigmay - Abs[sigma[t] - E0*epsi];
phi[t_, epsi_, dotepsi_] = (sigma[t] - E0 epsi) dotepsi;

events = {WhenEvent[g[t, epsi[t]] < phi[t, epsi[t], epsi'[t]], coef[t] -> 1], 
   WhenEvent[phi[t, epsi[t], epsi'[t]] < g[t, epsi[t]], 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, 100}, 
   DiscreteVariables -> coef, SolveDelayed -> True]

Plot[epsisol[t], {t, 0, 15}]

Answer 2

Ein Trick, um das volle Ergebnis zu erzielen.

sigma[t_] := Sin[t];
sigmay = 0.5;
E0 = 1;
tmax = Pi;
g[t_?NumericQ, epsi_] := sigmay - Abs[sigma[t] - E0*epsi]
phi[t_?NumericQ, epsi_, dotepsi_] := (sigma[t] - E0*epsi)*dotepsi
tmax = Pi;
tmin = 0;
epsisolant = sigma[tmin];
GR = {};

While[tmax < 100, 
  epsisol = NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, epsi[tmin] == epsisolant}, epsi, {t, tmin, tmax}, Method -> {"EquationSimplification" -> "Residual"}];
  AppendTo[GR, Plot[epsisol[t], {t, tmin, tmax}]];
  epsisolant = epsisol[tmax];
  tmin = tmax;
  tmax += Pi/4
]

Show[GR, PlotRange -> All]

Answer 3

Ihre Eingabe ist meines Wissens korrekt darauf angewendet. Gut gemacht.

Dies ist jedoch ein diskretisierter Versuch, das Problem zu lösen.

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, 10^-13, 100}]

Die zweite Nachricht öffnet eine Seite ndsolve :: ndcf mit der direkten Einladung, sich an den technischen Support von Wolfram Inc. zu wenden.

Ich fand heraus, dass die Domain mit einem schnellen Sprung von der Startzeit bei etwas mehr als 10 ^ -13 abhängt, zum Beispiel ein Viertel, in dem ich Ihre Ergebnisse reproduziere und um diese herum wieder weniger. Möglicherweise entspricht dies der Domainlänge 4.71. Das kann sogar bei höheren Startzeiten wie 0,0001 oder so wieder erreicht werden.

Meine Ausgabe ist:

Plot[epsisol[t], {t, 0.005, 4.71}, PlotRange -> Full]

Von da an stimme ich [ @cesareo ] 5 zu, dass es quasi periodisch verzögert weitergehen könnte. Dies könnte nicht nur in der Startzeit, sondern auch in der Zeit bereits chaotisch sein. Der Anstieg und Abfall kann charakteristisch sein. Irgendwie ähnelt das einem Sägezahn. Daher und weil die Switch-Funktion dies vorschlägt, stelle ich mir die Lösungsidee vor, dass dies mit Fourier- oder Laplace-Methoden für mehr Domain gelöst werden könnte. Dies funktioniert nur in Annäherung.

Aber die merkwürdige Idee hat meine Pläne geändert: die Domain willkürlich verkleinern:

epsisol = 
 NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, 
   epsi[0] == 0}, epsi, {t, 10^-13, 10}]

Plot[epsisol[t], {t, 0.005, 10}, PlotRange -> Full]

Hoffentlich hilft das. Dies geschieht mit V12.0.0 auf dem iMac Catalina.

Dies kann bis zu 10.99639 gelöst werden, wenn das Method -> {"EquationSimplification" -> "Residual"}verwendet wird. Die Nachricht bleibt : ndcf. Der wiederholte Konvergenztest akzeptiert nicht die schnelle Stagnation des Wachstums der Lösung bei -0,5. Es reicht jedoch für den gesamten Zeitraum des Diagramms aus. Vielleicht die Behandlung als differentiell-algebraische Gleichung.

Es scheint, dass eine bessere Übereinstimmung zwischen Sigmay und Sigma einen längeren Bereich in den Funktionen für differenziell-algebraische Standardmethoden bietet. Vielleicht ist dies andererseits eine Frage, die darauf ausgelegt ist, die Anpassungsfähigkeit des wiederholten Konvergenztests nicht zu erfüllen.

Ich habe einen Vergleich zwischen meiner und meiner Lösung von @xzczd gemacht .

Plot[{epsisol[t], epsisolu[t], 
  1.25 10^7 (epsisol[t] - epsisolu[t])}, {t, 0.00001, 10.99}, 
 PlotRange -> Full]

Obwohl beide Lösungen auf den ersten Blick sehr ähnlich aussehen, sind sie unterschiedlich.

Plot[{epsisol[t] - .5, epsisolu[t] - .5}, {t, 1.8, 3.14}, 
 PlotRange -> Full, PlotLegends -> "Expressions"]

Meins bleibt ein wenig, eine Größenordnung weiter von der Grenze von 0,5 entfernt. Dies ist für den negativen Rand noch größer und für das zweite konstante Intervall größer. Dann scheitert meine Lösung. Der selbst sehr kleine Fehler schwingt auf und am Ende schlägt der Test fehl.

Mit InterpolationOrder->Allder Schwingung werden viel kleiner und repetitiver:

Die Domain ist aber nicht größer.

Zu jedem Accuracy 9,10,11,... Es gibt ein Intervall nahe Null, für das die Integration erfolgreich ist.

epsisol = NDSolveValue [{Min [gi [t, epsi [t]], phi [t, epsi [t], epsi '[t]] == 0, epsi [0] == 0}, epsi, {t , 10 ^ -10.1295, 11}, Methode -> {"EquationSimplification" -> "Residual"}, InterpolationOrder -> All, AccuracyGoal -> 10]

Plot[{epsisol[t], epsisolu[t], 
  1.25 10^7 (epsisol[t] - epsisolu[t])}, {t, 0.00001, 7.85}, 
 PlotRange -> Full]

Für Accuracy 11 hat die Domäne ein sehr großes Intervall, für das meine Lösung der Referenzlösung viel näher kommt und die Schwingung gezähmt wird. Bei -0.5mir ist besser als bei dem Konkurrenten. Die Schwingungen bleiben aber immer noch in der Größenordnung von 10 ^ -7.

Schnell und schmutzig wie Mathematica heutzutage. Das Verhalten ist ein klarer Hinweis darauf, dass Mathematica StiffnessSwitching intern für die Berechnung der Lösung verwendet.

ListLinePlot@
 Quiet@Table[(epsisol = 
      NDSolveValue[{Min[gi[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 
         0, epsi[0] == 0}, epsi, {t, 10^expon, 11}, 
       Method -> {"EquationSimplification" -> "Residual"}, 
       InterpolationOrder -> All, AccuracyGoal -> 13])[[1, 1, 
     2]], {expon, -5, -16, -.01}]

Es gibt viele mögliche Startwerte für Genauigkeit 12. Das Ergebnis ist immer noch ein wildes Umschalten der Steifheit, aber die Genauigkeit nimmt stark zu.