Resolvendo um problema complementar não linear (plasticidade)

Question

Eu gostaria de resolver o seguinte: dado $t\mapsto\sigma(t)$ e $E>0$, $\sigma_y>0$, encontrar $\epsilon$ de tal modo 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.$$

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

Exemplo:

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

Tente com 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 *)

---

Contexto e detalhes

Estou interessado em plasticidade, em particular o cálculo do deslocamento $\varepsilon$ induzido por uma força aplicada $\sigma$ a um sistema composto por um controle deslizante e uma mola em paralelo, como:

Este problema é abordado na resolução de um sistema de equação diferencial-algébrica descontínua para o comportamento de plasticidade ($C_2$ é $H$ e $C_1$ é $\infty$) --- Acredito que haja um pequeno erro nas equações, mas ainda funciona após a correção.

Porém, com as duas respostas, não consegui adaptar o código para que funcionasse com os dois dispositivos em série. Além disso, eu queria derivar a equação da teoria da plasticidade "padrão":

• um critério de plasticidade $f =|\sigma -X| - \sigma_y \leq 0$ ($g=-f$ acima, e $X=E \epsilon$ é o estresse na primavera);
• a positividade da dissipação: $\phi = \epsilon' (\sigma_y -E\epsilon)\geq 0$
• uma condição de ortogonalidade que implica que a energia é dissipada se houver plasticidade ($f=0$): $f\times \phi = 0$.

Muitas vezes isso é escrito completamente: $$ 0\leq (-f) \perp \phi \geq 0$$

Esses tipos de formulação também são encontrados na dinâmica de contato intermitente: a força de reação é sempre não negativa, e pode ser diferente de zero apenas se houver contato, ou seja, quando a distância for zero. Por outro lado, se a distância não for negativa, a reação para só pode ser zero.

Esses problemas são numericamente desafiadores, embora existam métodos numéricos dedicados. Mesmo a formulação envolvendo o derivado$\epsilon'$ é vacilante, porque $\epsilon$ não é diferenciável em todos os lugares (assim como a velocidade de uma bola quicando não é definida no momento do impacto).

De qualquer forma, WhenEventfunciona muito bem para bolas quicando com poucos contatos, então eu esperava WhenEventser eficiente aqui.

Answer 1

Acho que o comportamento do WhenEventencontrado pelo OP é um bug. Enfim, aqui está uma WhenEventsolução baseada em trabalho :

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

Um truque para obter o resultado completo.

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

A sua contribuição é, pelo meu conhecimento, aplicada a ela corretamente. Bem feito.

Mas esta é uma tentativa discretizada de resolver o problema.

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

A segunda mensagem abre uma página ndsolve :: ndcf com o convite direto para entrar em contato com o suporte técnico da Wolfram Inc.

Descobri que o domínio depende com um salto rápido do tempo de início um pouco mais de 10 ^ -13 por exemplo, um quarto eu reproduzo Seus resultados e menos novamente. Pode ser que seja uma correspondência para o comprimento do domínio 4,71. Isso pode até ser obtido novamente em tempos de início mais altos, como 0,0001 ou mais.

Minha saída é:

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

A partir daí, eu concordo com [ @cesareo ] 5 , pode continuar atrasado quase periódico. Isso já pode ser caótico não apenas no horário de início, mas também no período. A ascensão e queda podem ser características. De alguma forma, isso é semelhante a um dente de serra. Portanto, e porque a função switch sugere isso, faço a ideia de solução de que isso pode ser resolvido com as metodologias de Fourier ou Laplace para mais domínio. Isso só funcionará em aproximação.

Mas a curiosa ideia mudou meus planos: tornar o domínio menor arbitrariamente:

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]

Espero que ajude. Isso é feito com a V12.0.0 no iMac Catalina.

Isso pode ser resolvido em até 10.99639 se o Method -> {"EquationSimplification" -> "Residual"}for usado. A mensagem permanece: ndcf. O teste de convergência repetido não aceita a estagnação rápida do crescimento da solução em -0,5. Mas é suficiente para todo o período do gráfico. Talvez o tratamento como uma equação diferencial-algébrica.

Parece que uma correspondência melhor para sigmay e sigma dá um domínio mais longo nas capacidades dos métodos algébricos diferenciais de prateleira. Talvez esta seja, por outro lado, uma questão projetada para a falha de adaptabilidade do teste de convergência repetido.

Fiz uma comparação entre as duas soluções, a minha e a da @xzczd .

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

Apesar de ambas as soluções parecerem à primeira vista muito semelhantes, são diferentes.

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

O meu fica um pouco, uma ordem de magnitude além do limite de 0,5. Isso é ainda maior para a borda negativa e maior no segundo intervalo constante. Então minha solução vai ao fracasso. Mesmo o erro muito pequeno oscila para cima e no final o teste falha.

Com InterpolationOrder->Alla oscilação fica muito menor e mais repetitiva:

Mas o domínio não é maior.

Para cada um Accuracy 9,10,11,... há um intervalo próximo a zero para o qual a integração é bem-sucedida.

epsisol = NDSolveValue [{Min [gi [t, epsi [t]], phi [t, epsi [t], epsi '[t]]] == 0, epsi [0] == 0}, epsi, {t , 10 ^ -10,1295, 11}, Método -> {"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]

Para Precisão 11, o domínio tem um intervalo muito grande para o qual minha solução fica muito mais próxima da solução de referência e a oscilação é domada. No -0.5meu é melhor do que no concorrente. Mas as oscilações permanecem na ordem de 10 ^ -7.

Rápido e sujo como os built-ins do Mathematica são hoje em dia. O comportamento é uma dica clara de que o Mathematica usa o StiffnessSwitching internamente para o cálculo da solução.

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

Existem muitos valores iniciais possíveis para a Precisão 12. O resultado ainda é a rigidez alternando descontroladamente, mas a precisão aumenta fortemente.