Résolution d'un problème complémentaire non linéaire (plasticité)

Question

Je voudrais résoudre ce qui suit: étant donné $t\mapsto\sigma(t)$ et $E>0$, $\sigma_y>0$, trouver $\epsilon$ tel que:

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

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

Exemple:

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

Tentative avec 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 *)

---

Contexte et détails

Je m'intéresse à la plasticité, en particulier au calcul du déplacement $\varepsilon$ induit par une force appliquée $\sigma$ à un système composé d'un curseur et d'un ressort en parallèle, tel que:

Ce problème est traité dans Résolution d'un système d'équations différentielles-algébriques discontinues pour le comportement de plasticité ($C_2$ est $H$ et $C_1$ est $\infty$) --- Je crois qu'il y a une légère erreur dans les équations mais cela fonctionne toujours après correction.

Cependant, avec les deux réponses, je n'ai pas réussi à adapter le code pour qu'il fonctionne avec deux appareils de ce type en série. Aussi, je voulais dériver l'équation de la théorie de plasticité "standard":

• un critère de plasticité $f =|\sigma -X| - \sigma_y \leq 0$ ($g=-f$ ci-dessus, et $X=E \epsilon$ est la contrainte au printemps);
• la positivité de la dissipation: $\phi = \epsilon' (\sigma_y -E\epsilon)\geq 0$
• une condition d'orthogonalité impliquant que l'énergie est dissipée ssil y a plasticité ($f=0$): $f\times \phi = 0$.

Ceci est souvent écrit complètement: $$ 0\leq (-f) \perp \phi \geq 0$$

De tels types de formulation se retrouvent également dans les dynamiques de contact intermittentes: la force de réaction est toujours non négative, et ne peut être non nulle que s'il y a contact, c'est-à-dire lorsque la distance est nulle. A l'inverse, si la distance n'est pas négative, la réaction pour ne peut être que nulle.

Ces problèmes sont numériquement difficiles, même s'il existe des méthodes numériques dédiées. Même la formulation impliquant le dérivé$\epsilon'$ est bancal, car $\epsilon$ n'est pas différentiable partout (tout comme la vitesse d'une balle rebondissante n'est pas définie au moment de l'impact).

Quoi qu'il en soit, WhenEventfonctionne très bien pour faire rebondir des balles avec peu de contacts, donc je me serais attendu WhenEventà être efficace ici.

Answer 1

Je pense que le comportement WhenEventrencontré par OP est un bug. Quoi qu'il en soit, voici une WhenEventsolution fonctionnelle:

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

Une astuce pour obtenir le résultat complet.

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

Votre entrée est à ma connaissance appliquée correctement. Bien joué.

Mais il s'agit d'une tentative discrétisée de résoudre le problème.

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}]

Le deuxième message ouvre une page ndsolve :: ndcf avec l'invitation directe à contacter le support technique de Wolfram Inc.

J'ai trouvé que le domaine dépend d'un saut rapide sur l'heure de départ à un peu plus de 10 ^ -13 par exemple un quart je reproduis vos résultats et autour de cela moins encore. Il se peut que ce soit une correspondance pour la longueur de domaine 4.71. Cela peut même être obtenu à des heures de départ plus élevées que 0,0001 ou plus.

Ma sortie est:

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

A partir de là, je suis d'accord avec [ @cesareo ] 5, cela pourrait être retardé quasi-périodique. Cela pourrait déjà être chaotique non seulement dans le temps de départ mais dans la période. La montée et la descente peuvent être caractéristiques. D'une manière ou d'une autre, cela ressemble à une dent de scie. Par conséquent et parce que la fonction de commutation le suggère, je fais l'idée de solution que cela pourrait être résolu avec les méthodologies de Fourier ou Laplace pour plus de domaine. Cela ne fonctionnera qu'en approximation.

Mais l'idée curieuse a changé mes plans: réduire arbitrairement le domaine:

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]

J'espère que ça t'as aidé. Ceci est fait avec V12.0.0 sur iMac Catalina.

Cela peut être résolu jusqu'à 10,99639 si le Method -> {"EquationSimplification" -> "Residual"}est utilisé. Le message reste: ndcf. Le test de convergence répété n'accepte pas la stagnation rapide de la croissance de la solution à -0,5. Mais cela suffit pour toute la période du graphique. Peut-être le traitement comme une équation différentielle-algébrique.

Il semble qu'une meilleure correspondance pour sigmay et sigma donne un domaine plus long dans les capacités des méthodes d'algèbre différentielle du marché. Peut-être s'agit-il d'autre part d'une question destinée à l'échec de l'adaptabilité du test de convergence répétée.

J'ai fait une comparaison entre la solution, la mienne et celle de @xzczd .

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

Malgré les deux solutions semblent à première vue très similaires, elles sont différentes.

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

Le mien reste un peu, un ordre d'une grandeur plus loin de la limite de 0,5. C'est encore plus grand pour la frontière négative et plus grand au deuxième intervalle constant. Ensuite, ma solution échoue. L'erreur même très petite oscille vers le haut et la fin est l'échec du test.

Avec InterpolationOrder->Alll'oscillation, devenez beaucoup plus petit et répétitif:

Mais le domaine n'est pas plus grand.

À chacun Accuracy 9,10,11,... il existe un intervalle proche de zéro pour lequel l'intégration est réussie.

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

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

Pour Accuracy 11, le domaine a un intervalle beaucoup plus grand pour lequel ma solution se rapproche beaucoup plus de la solution de référence et l'oscillation est apprivoisée. Chez -0.5moi, c'est mieux que celui du concurrent. Mais les oscillations restent toujours d'ordre 10 ^ -7.

Rapide et sale comme les intégrés Mathematica le sont ces jours-ci. Le comportement indique clairement que Mathematica utilise StiffnessSwitching en interne pour le calcul de la solution.

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}]

Il existe de nombreuses valeurs de départ possibles pour la précision 12. Le résultat est encore une variation brutale de la rigidité, mais la précision augmente fortement.