​

==Resumen==

​

El problema de la simulación de materiales no lineales sujetos a inestabilidad estructural implica en numerosos casos la divergencia de métodos implícitos que usan algoritmos de retorno para la simulación del comportamiento constitutivo. En este artículo se presenta una estrategia combinada entre algoritmos implícita y explícita para la solución de problemas mediante el método de elementos finitos que resuelve este inconveniente. No se considera la partición de la malla de elementos finitos, es decir, todo el dominio es resuelto implícita o explícitamente en un determinado instante. Se considera la aplicación a materiales hiperelásticos no lineales y a materiales elastoplásticos no lineales. Adicionalmente, la respuesta del sólido ante grandes deformaciones es formulada e integrada en la estrategia propuesta. Distintos ejemplos numéricos en problemas no lineales (geométricos y/o materiales) han demostrado la efectividad de la técnica.

​

==Abstract==

​

Simulation problems involving non-linear materials imply in numerous cases divergence of the implicit method which use return mapping algorithms for modelling of the nonlinear response. A switching implicit-explicit numerical technique in the context of Finite Element Methods is presented in this paper. Implicit/explicit mesh partitions are not considered whatsoever. Formulation for application to nonlinear hyperelastic materials and nonlinear elastic-plastic materials is provided. Furthermore, the response of the solid subjected to large deformations is presented and is embedded in the proposed technique. Numerical tests for nonlinear problems (geometric and/or material) showed the accurateness of the technique.

​

==Palabras clave==

​

Método de elementos finitos ; Inestabilidad ; Estructuras ; Implícito ; Explícito

​

==Keywords==

​

Finite element method ; Instability ; Structures ; Implicit ; Explicit

​

==1. Introducción==

​

Se ha propuesto un gran número de técnicas para la resolución de problemas de convergencia en el tratamiento y simulación de sólidos no lineales e inelásticos. El método de elementos finitos (MEF) y métodos derivados del MEF han demostrado ser, en general, procedimientos numéricos convenientes para la solución de problemas en elastodinámica y, en particular, para el análisis de deformaciones de sólidos no lineales. Sin embargo, el procedimiento de actualización de las tensiones en los puntos de cuadratura se lleva a cabo mediante algoritmos de retorno [[#bib0005|[1]]]  and [[#bib0010|[2]]] , los cuales conllevan con frecuencia problemas de estabilidad aparte de ser costosos. Puede consultarse en Kojić [[#bib0015|[3]]]  una excelente revisión de las técnicas para la resolución de materiales inelásticos.      

​

La idea propuesta en este trabajo consiste en empezar la ejecución con un algoritmo implícito, generalmente más rápido para este tipo de problemas y, si existe posibilidad de divergencia, cambiar a una ejecución explícita (condicionalmente estable). Si el origen de la divergencia desaparece, el retorno al algoritmo más rápido, en este caso el implícito, es una opción deseable en el caso de una ejecución costosa. En la literatura pueden encontrarse técnicas de este tipo aplicadas a diferencias finitas. Así, Ascher et al. [[#bib0020|[4]]]  proponen un método implícito-explícito para diferencias finitas centradas (aparentemente pionero) para la integración temporal de ecuaciones en derivadas parciales. Ruuth [[#bib0025|[5]]]  lo aplicó con éxito a la resolución de problemas de reacción-difusión. Más recientemente, Robinson [[#bib0030|[6]]]  ha propuesto un método implícito-explícito en diferencias finitas para ecuaciones parabólicas semilineales probando la existencia y unicidad de la solución. Idesman et al. [[#bib0035|[7]]]  aplican una discretización temporal implícita-explícita usando elementos finitos para el problema de propagación de ondas. Comblen et al. [[#bib0040|[8]]]  presentan un trabajo interesante para la simulación de modelos marinos usando una partición temporal aplicada a métodos Runge-Kutta.      

​

En este artículo se propone una técnica temporal implícita-explícita para problemas inelásticos con extensión a grandes deformaciones. No se incluyen en este trabajo las técnicas numéricas que consideran la partición de la malla de elementos finitos para resolución implícita o explícita, tal como las que presentan Hughes et al. [[#bib0045|[9]]]  o Farhat et al. [[#bib0050|[10]]] . La resolución es secuencial en el tiempo, es decir, el flujo de la solución puede volcarse a uno u otro resolvedor (implícito o explícito) según corresponda, dependiendo de las condiciones de estabilidad y de la estimación del error [[#bib0055|[11]]] . Otro tipo de problemas que podrían beneficiarse de la técnica propuesta, a pesar de que no se han intentado en este trabajo, pueden ser los problemas de contacto no suave que involucran esquinas agudas y una superficie cóncava, por ejemplo. Las condiciones no lineales de fricción sobre geometrías complejas han demostrado ser divergentes cuando se ha usado un resolvedor implícito; véanse por ejemplo los trabajos de Oñate y Flores [[#bib0060|[12]]]  o de Laursen y Simo [[#bib0065|[13]]] .      

​

En la siguiente sección de este artículo se introducen los puntos especiales de la curva de equilibrio que suelen crear divergencia en estructuras esbeltas cuando se usan modelos no-lineales en conjunción con un método implícito y para los cuales la técnica propuesta es especialmente efectiva. En segundo lugar, se enseña el modelo elastoplástico implementado en el MEF codificado con atención a la actualización de las tensiones usando algoritmos de retorno. A continuación, se muestra el problema incremental resuelto implícita o explícitamente incluyendo las condiciones de transferencia de flujo. Finalmente, se presentan test numéricos mostrando la robustez de la técnica propuesta. Estos incluyen ''snap-through''  en un arco esbelto, arco circular elastoplástico y grandes deformaciones en la membrana de Cook.      

​

==2. Inestabilidad estructural==

​

Desde un punto de vista matemático, la inestabilidad estructural se caracteriza por la aparición de puntos especiales en la curva de equilibrio. Estos puntos especiales pueden ser agrupados en 2 [[#bib0070|[14]]] :

* Puntos de retorno, los cuales provocan en la estructura el abandono del equilibrio debido a fallos estructurales del material. Puede ocurrir que la estructura alcance un nuevo estado de equilibrio o que, por el contrario, conlleve el colapso total de la estructura.

* Puntos críticos (límites), los cuales suelen afectar al método de solución. El equilibrio entre puntos límites no ofrece ningún problema a los métodos usados. Sin embargo, pueden aparecer inestabilidades numéricas en los puntos límites, haciendo diverger al método. Los modos más característicos en los que aparecen son:

* ''Snap-through:''  la curva de equilibrio es no lineal, a partir del primer punto límite se produce un reblandecimiento hasta el segundo punto límite. Tiene lugar un endurecimiento a medida que la deformación aumenta a partir de ese segundo punto límite. Este tipo de comportamiento puede ser observado en arcos esbeltos cargados en su punto medio.                                      

* ''Snap-back:''  en este caso, la curva de equilibrio es dirigida hacia atrás (incremento de deformación negativo) después de alcanzar el primer punto límite.                              

​

Los problemas de ''snap-through''  se han resuelto con la técnica propuesta, ya que este tipo de problemas no converge con Newton-Raphson a menos que se usen procedimientos ''arc-length'' . El lector interesado en este tipo de métodos, y sus inconvenientes, puede consultar las referencias  [[#bib0075|[15]]] , [[#bib0080|[16]]] , [[#bib0085|[17]]]  and [[#bib0090|[18]]] . Matías y Oñate [[#bib0095|[19]]]  presentan otro enfoque mediante el método de desplazamientos críticos, en el cual se predice un camino crítico aproximado imponiendo la condición de singularidad mediante una expresión aproximada de la matriz de rigidez tangente en el punto crítico.      

​

==3. Modelo constitutivo elastoplástico==

​

En esta sección se expone brevemente la base de los modelos elastoplásticos implementados en el programa para esta técnica. Se muestra en la sección [[#sec0035|4]]  un ejemplo detallado para el caso específico del modelo de Von-Mises. Los elementos básicos de un modelo constitutivo elastoplástico son:

* Descomposición de la deformación en elástica y plás-tica.

* Evolución elástica.

* Criterio de límite elástico, el cual es representado por una superficie en el espacio de tensiones.

* Evolución de la deformación plástica: regla de flujo plástico.

* Evolución de la superficie de límite elástico: ley de endurecimiento.

​

===3.1. Límite elástico===

​

Esta superficie se define generalmente en el espacio de tensiones delimitando el dominio elástico. La expresión que la representa, ''Y'' , toma un valor negativo para deformaciones elásticas y valor nulo en el momento en el que tiene lugar la deformación plástica. El dominio elástico ''E''  se define como:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>E=\lbrace \boldsymbol{\sigma }\parallel Y(\boldsymbol{\sigma },\boldsymbol{\alpha })\quad <0\rbrace </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)

|}

​

donde '''''σ'''''  es el tensor de tensiones y '''''α'''''  es el tensor que contiene las variables de estado internas del material. La superficie de límite elástico ''P''  es:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>P=\lbrace \boldsymbol{\sigma }\parallel Y(\boldsymbol{\sigma },\boldsymbol{\alpha })=</math><math>0\rbrace </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)

|}

​

===3.2. Regla de flujo plástico y ley de endurecimiento===

​

El modelador decide la evolución de la deformación plástica y de las variables de estado. Así, los trabajos debidos al endurecimiento o a la deformación plástica acumulada son generalmente aceptados como variables internas de estado. En el caso de la existencia de daño, se necesita una variable entre cero y la unidad que degrade la rigidez del material. La regla de flujo plástico y la ley de endurecimiento pueden ser postuladas como sigue:

​

<span id='eq0015'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\dot{\epsilon }}}^p=\dot{\gamma }\quad \boldsymbol{\mbox{N}}(\boldsymbol{\sigma },\boldsymbol{\alpha })=</math><math>\dot{\gamma }\quad \frac{\partial \Psi }{\partial \boldsymbol{\sigma }}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)

|}

​

donde '''N''' ('''''σ''''' , '''''α''''' ) es el vector de flujo. La ley de endurecimiento es:

​

<span id='eq0020'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\dot{\alpha }}=\dot{\gamma }\quad \boldsymbol{\mbox{H}}(\boldsymbol{\sigma },\boldsymbol{\alpha })</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)

|}

​

donde '''H''' ('''''σ''''' , '''''α''''' ) define las variables de endurecimiento y es llamado módulo de endurecimiento generalizado [[#bib0100|[20]]] . '''N'''  y '''H'''  pueden derivarse desde el potencial de flujo. En el caso de que este potencial sea la superficie de límite elástico, el flujo se denomina asociativo. Los metales son modelados, generalmente, de flujo asociativos.      

​

===3.3. Criterio de carga/descarga===

​

La evolución de flujo plástico, ecuación [[#eq0015|(3)]] , y de endurecimiento, ecuación [[#eq0020|(4)]] , se complementa con el criterio de carga/descarga, representado por las condiciones de Kuhn-Tucker:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Phi \quad \leq \quad 0\quad \dot{\gamma }\geq \quad 0\quad \Phi \quad \dot{\gamma }=</math><math>0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)

|}

​

<span id='sec0035'></span>

==4. Modelo constitutivo de Von-Mises==

​

===4.1. Criterio de límite elástico de Von-Mises===

​

En esta sección se recuerda brevemente el modelo de Von-Mises para el cual se implementó uno de los algoritmos de retorno. Este criterio postula la plasticidad del material cuando el segundo invariante del tensor de tensiones desviador ''J''<sub>2</sub> ('''s''' ) alcanza un valor crítico ''σ''<sup>''y''</sup> . El tensor desviador '''s'''  es función de las variables de estado. La presión hidrostática no influye en el criterio tal y como ocurre con el criterio de Tresca. En un estado de cortadura pura, es decir, círculo de Mohr centrado en el origen, ''σ''<sub>1</sub>  = − ''σ''<sub>2</sub>  > 0, ''σ''<sub>3</sub>  = 0 y <math display="inline">\sqrt{J_2(\boldsymbol{\mbox{s}})}=</math><math>{\tau }_{max}={\sigma }_1</math>  la función toma la forma:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Phi (\boldsymbol{\sigma })=\sqrt{J_2(\boldsymbol{\mbox{s}})}-</math><math>{\tau }^y</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)

|}

​

En el caso de tensiones unidireccionales, el criterio se postula:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Phi (\boldsymbol{\sigma })=\sqrt{3\quad J_2(\boldsymbol{\mbox{s}})}-</math><math>{\sigma }^y</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)

|}

​

donde ''σ''<sup>''y''</sup>  es el valor del límite elástico obtenido en el ensayo de tracción. <math display="inline">\sqrt{3\quad J_2(\boldsymbol{\mbox{s}})}</math>  es conocido como la tensión efectiva de Von-Mises o tensión equivalente. Véase que:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\sigma }^y=\sqrt{3}\quad {\tau }^y</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

La superficie de Von-Mises en el caso general tridimensional es configurada como un cilindro infinito con la linea hidrostática como eje de dicho cilindro.

​

==5. Algoritmo de integración para el modelo de Von-Mises isótropo con endurecimiento==

​

En un intérvalo genérico de tiempo [''t''<sub>''n''</sub> , ''t''<sub>''n'' +1                            </sub> ], el incremento de la tasa de deformaciones es dado por la ecuación [[#eq0045|(8)]] . Las variables de estado consideradas son la deformación elástica <math display="inline">{\boldsymbol{\epsilon }}_n^e</math>  y la deformación plástica acumulada <math display="inline">{\boldsymbol{\overline{\epsilon }}}_n^p</math>  al principio del intervalo de tiempo [''t''<sub>''n''</sub> , ''t''<sub>''n'' +1                            </sub> ]. La predicción del estado elástico de deformaciones viene dada por la ecuación [[#eq0050|(9)]] . La ecuación [[#eq0055|(10)]]  recoge que no se incremente la deformación plástica:

​

<span id='eq0045'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta \boldsymbol{\epsilon }={\boldsymbol{\epsilon }}_{n+1}-</math><math>{\boldsymbol{\epsilon }}_n</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)

|}

​

<span id='eq0050'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\epsilon }}_{n+1}^{etrial}={\boldsymbol{\epsilon }}_n^e+</math><math>\Delta \boldsymbol{\epsilon }</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)

|}

​

<span id='eq0055'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\overline{\epsilon }}}_{n+1}^{ptrial}=</math><math>{\boldsymbol{\overline{\epsilon }}}_n^p</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)

|}

​

El tensor de tensiones predecido asumiendo una evolución plástica viene dado por la ecuación [[#eq0060|(11)]] . La tensión correspondiente al límite elástico depende de la deformación plástica acumulada en el instante ''t''<sub>''n''</sub> , ecuación [[#eq0065|(12)]] :

​

<span id='eq0060'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\sigma }}_{n+1}^{trial}=\boldsymbol{\mbox{D}}^e:{\boldsymbol{\epsilon }}_{n+1}^{etrial}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)

|}

​

<span id='eq0065'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\sigma }_{yn+1}^{trial}={\sigma }_y({\boldsymbol{\overline{\epsilon }}}_n^p)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)

|}

​

Si <math display="inline">{\boldsymbol{\sigma }}_{n+1}^{trial}</math>  yace en el interior del dominio elástico, ecuación [[#eq0070|(13)]] , la evolución es puramente elástica en el intervalo de tiempo [''t''<sub>''n''</sub> , ''t''<sub>''n'' +1                            </sub> ] y la predicción corresponde a la solución para ese incremento:

​

<span id='eq0070'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Phi ({\boldsymbol{\sigma }}_{n+1}^{trial},{\sigma }_{yn})\leq 0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)

|}

​

La actualización es simplemente llevada a cabo como sigue:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\epsilon }}_{n+1}^e={\boldsymbol{\epsilon }}_{n+1}^{etrial}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\overline{\epsilon }}}_{n+1}^p={\boldsymbol{\overline{\epsilon }}}_{n+1}^{ptrial}=</math><math>{\boldsymbol{\overline{\epsilon }}}_n^p</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\sigma }}_{n+1}={\boldsymbol{\sigma }}_{n+1}^{trial}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\sigma }_{yn+1}={\sigma }_{yn+1}^{trial}={\sigma }_{yn}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)

|}

​

En caso contrario, la evolución es elastoplástica y el estado predecido yace fuera del contorno del dominio elástico. En este caso, el mapeo de retorno (véase [[#fig0005|fig. 1]] ) es necesario, como sigue:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\epsilon }}_{n+1}^e={\boldsymbol{\epsilon }}_{n+1}^{etrial}-</math><math>\Delta \gamma \sqrt{\frac{2}{3}}\frac{\boldsymbol{\mbox{s}}_{n+1}}{\parallel \boldsymbol{\mbox{s}}_{n+1}\parallel }</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\overline{\epsilon }}}_{n+1}^p={\boldsymbol{\overline{\epsilon }}}_n^p+</math><math>\Delta \gamma </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>0=\sqrt{3\quad J_2(\boldsymbol{\mbox{s}}_{n+1})}-{\sigma }_y({\boldsymbol{\overline{\epsilon }}}_{n+1}^p)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)

|}

​

<span id='fig0005'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_721884677-1-s2.0-S0213131513000217-gr1.jpg|center|356px|Pseudo integración para caracterización del endurecimiento: predicción elástica ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 1.

​

Pseudo integración para caracterización del endurecimiento: predicción elástica – corrección plástica.

​

</span>

|}

​

Este grupo de ecuaciones algebraicas tiene que ser resuelto para <math display="inline">{\boldsymbol{\epsilon }}_{n+1}^e</math> , <math display="inline">{\boldsymbol{\overline{\epsilon }}}_{n+1}^p</math>  y Δ''γ'' . '''s'''<sub>''n'' +1                            </sub>  es el tensor desviador en ''t''<sub>''n'' +1                            </sub>  ecuación [[#eq0110|(21)]] :

​

<span id='eq0110'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{s}}_{n+1}=2Gdev[{\boldsymbol{\epsilon }}_{n+1}^e]</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)

|}

​

Este sistema puede ser reducido a una sola ecuación, según Souza [[#bib0105|[21]]] , teniendo como variable al multiplicador plástico Δ''γ'' . Esta contracción del sistema hace el procedimiento menos costoso.      

​

==6. Algoritmo implícito==

​

En este caso se usa el método de Newton-Raphson para resolver el sistema producido por la forma débil de las ecuaciones de momento. El algoritmo de retorno es invocado en cada incremento del proceso general para la actualización de tensiones como se muestra en el cuadro I. El algoritmo implícito se basa en una pseudo discretización temporal siguiendo el trabajo de Simo y Hughes [[#bib0010|[2]]]  considerando, así, la transición de deformación entre 2 puntos de tiempo. Concretamente, se ha implementado el método de Euler hacia atrás en conjunción con el método iterativo de Newton-Raphson [[#bib0100|[20]]] . Así, si en un incremento de tiempo [''t''<sub>''n''</sub> , ''t''<sub>''n'' +1                            </sub> ] se da un conjunto de variables internas ''α''<sub>''n''</sub>  en el instante ''t''<sub>''n''</sub> , el tensor de deformaciones '''''ε''''' (''t''<sub>''n'' +1                            </sub> ) debe determinar las tensiones '''''σ''''' (''t''<sub>''n'' +1                            </sub> ) y las variables internas solo a través del algoritmo de integración. Véase Simo y Hughes [[#bib0010|[2]]]  para detalles de este tipo de algoritmos, es decir:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\sigma }(t_{n+1})=\boldsymbol{\overset{\mbox{ˆ}}{\sigma }}({\boldsymbol{\alpha }}_n,{\boldsymbol{\epsilon }}_{n+1})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\alpha }(t_{n+1})=\boldsymbol{\overset{\mbox{ˆ}}{\alpha }}({\boldsymbol{\alpha }}_n,{\boldsymbol{\epsilon }}_{n+1})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)

|}

​

Después de la discretización mediante elementos finitos, el problema se centra en encontrar los desplazamientos nodales '''''u'''''<sub>''n'' +1                            </sub>  en el instante ''t''<sub>''n'' +1                            </sub> , de modo que se satisfaga el sistema de ecuaciones incremental no lineal [[#eq0125|(24)]] :

​

<span id='eq0125'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{R}(\boldsymbol{u}_{n+1})=\boldsymbol{\mbox{f}}^{int}(\boldsymbol{u}_{n+1})-</math><math>\boldsymbol{\mbox{f}}_{n+1}^{ext}=0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)

|}

​

donde las fuerzas nodales internas '''f'''<sup>''int''</sup> ('''''u'''''<sub>''n'' +1                            </sub> ) y externas <math display="inline">\boldsymbol{\mbox{f}}_{n+1}^{ext}</math>  se obtienen como sigue:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{f}}^{int}(\boldsymbol{u}_{n+1})=</math><math>{\bigwedge}_{e=1}^{nelem}\left\{{\int }_{{\Omega }^{\left(e\right)}}\boldsymbol{\mbox{B}}^T\boldsymbol{\overset{\mbox{ˆ}}{\sigma }}({\boldsymbol{\alpha }}_n,\boldsymbol{\epsilon }(\boldsymbol{u}_{n+1}))\quad dv\right\}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{f}}_{n+1}^{ext}={\bigwedge}_{e=1}^{nelem}\left\{{\int }_{{\Omega }^{\left(e\right)}}\boldsymbol{\mbox{N}}^T\boldsymbol{b}_{n+1}dv+\right. </math><math>\left. {\int }_{\partial {\Omega }^{\left(e\right)}}\boldsymbol{\mbox{N}}^T\quad \boldsymbol{q}_{n+1}ds\right\}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)

|}

​

donde '''N'''  contiene las ''N'' (''ξ'' , ''η'' ) funciones de formas bilineales, '''b'''<sub>''n'' +1                            </sub>  son las fuerzas volumétricas, '''q'''<sub>''n'' +1                            </sub>  las fuerzas actuando sobre el contorno del sólido y '''B'''  el operador de deformaciones el cual tiene el siguiente formato en tensión/deformación plana para un elemento genérico (''e'' ) (el primer índice denota el número de nodo y la coma derivada):

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{B}}=\left[\begin{array}{ccccccc}

N_{1,1}^{\left(e\right)} & 0 & N_{2,1}^{\left(e\right)} & 0 & \ldots  & N_{n_{node},1}^{\left(e\right)} & 0\\

0 & N_{1,2}^{\left(e\right)} & 0 & N_{2,2}^{\left(e\right)} & \ldots  & 0 & N_{n_{node},2}^{\left(e\right)}\\

N_{1,2}^{\left(e\right)} & N_{1,1}^{\left(e\right)} & N_{2,2}^{\left(e\right)} & N_{2,1}^{\left(e\right)} & \ldots  & N_{n_{node},2}^{\left(e\right)} & N_{n_{node},1}^{\left(e\right)}\\

 &  &  &  &  &  & 

\end{array}\right]</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

Para continuar con el procedimiento numérico, la ecuación [[#eq0125|(24)]]  necesita ser linealizada. Pueden encontrarse en Souza et al. [[#bib0100|[20]]]  los detalles de esta operación.      

​

==7. Solución del problema incremental==

​

Como se ha mencionado, el método de Newton-Raphson se usa como parte de la estrategia implícita [véase la ecuación [[#eq0125|(24)]] ]. En el caso de materiales no lineales, en cada iteración hay que resolver la versión linealizada [ecuación [[#eq0125|(24)]] ] para el cálculo de desplazamientos nodales ''δ'''''''u'''''<sup>(''k'' )                            </sup> :

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{K}}_T\delta \boldsymbol{u}^{\left(k\right)}=</math><math>-\boldsymbol{R}^{\left(k-1\right)}(\boldsymbol{u}_{n+1})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)

|}

​

donde '''K'''<sub>''T''</sub>  es la matriz tangente de rigidez dada por:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{K}}_T=\frac{\partial \boldsymbol{R}}{\partial \quad \boldsymbol{u}_{n+1}}\vert _{\boldsymbol{u}_{n+1}^{\left(k-1\right)}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)

|}

​

la cual es obtenida mediante ensamblaje de las matrices tangentes de rigidez de cada elemento:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{k}}_T^{\left(e\right)}={\int }_{{\Omega }^{\left(e\right)}}\boldsymbol{\mbox{B}}^T\quad \overset{\mbox{ˆ}}{\boldsymbol{\mbox{D}}}\boldsymbol{\mbox{B}}dv</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)

|}

​

​

[[Image:draft_Content_721884677-1-s2.0-S0213131513000217-fx1.jpg|center|383px|Image for unlabelled figure]]

​

donde <math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{\mbox{D}}}</math>  es la matriz consistente tangente de rigidez [[#bib0100|[20]]] :

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\overset{\mbox{ˆ}}{\boldsymbol{\mbox{D}}}=\frac{\partial \boldsymbol{\overset{\mbox{ˆ}}{\sigma }}}{\partial \boldsymbol{{\epsilon }_{n+1}}}\quad \vert _{{\boldsymbol{\epsilon }}_{n+1}^{\left(k-1\right)}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)

|}

​

La estrategia global se muestra con detalle en el cuadro I. Obsérvese que, si la versión implícita diverge, el flujo de la ejecución se desvía a la versión explícita en el paso 7. Las condiciones de transición entre flujo implícito y explícito se explican en la sección [[#sec0090|9]] .      

​

===7.1. Caso de materiales sujetos a grandes deformaciones===

​

En el caso de elasticidad lineal, la respuesta o el comportamiento del material es independiente del camino de carga seguido, por tanto, no hay necesidad de usar un algoritmo de retorno para la actualización de su estado interno. Sin embargo, esto no es posible para el caso de materiales no lineales cuyo comportamiento se ve afectado por la historia de carga ejercida. Adicionalmente, el caso de grandes deformaciones es generalmente dependiente del camino y, por tanto, es necesario el uso de un algoritmo de retorno en cada iteración para alcanzar la solución.

​

Las tensiones se representan en este caso por <math display="inline">{\boldsymbol{\sigma }}_{n+1}=</math><math>\boldsymbol{\overset{\mbox{ˆ}}{\sigma }}({\boldsymbol{\alpha }}_n,\boldsymbol{\mbox{F}}_{n+1})</math> , donde la parte derecha corresponde a la tensión obtenida mediante el algoritmo de retorno. '''F'''<sub>''n'' +1                            </sub>  es el gradiente de deformaciones calculado al final del intervalo [''t''<sub>''n''</sub> , ''t''<sub>''n'' +1                            </sub> ]. Ahora los vectores que contienen las fuerzas nodales (internas y externas) están basados en la configuración deformada:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{f}}^{int}(\boldsymbol{u}_{n+1})=</math><math></math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\bigwedge}_{e=1}^{nelem}\left\{{\int }_{{\varphi }_{n+1}({\Omega }^{\left(e\right)})}\boldsymbol{\mbox{B}}^T\boldsymbol{\sigma }({\boldsymbol{\alpha }}_n,\boldsymbol{\mbox{F}}(\boldsymbol{u}_{n+1}))dv\right\}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{f}}_{n+1}^{ext}={\bigwedge}_{e=1}^{nelem}{\int }_{{\varphi }_{n+1}({\Omega }^{\left(e\right)})}\boldsymbol{\mbox{N}}^T\boldsymbol{b}_{n+1}dv+</math><math></math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>+{\int }_{\partial {\varphi }_{n+1}({\Omega }^{\left(e\right)})}\boldsymbol{\mbox{N}}^T\boldsymbol{q}_{n+1}ds</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)

|}

​

donde ''φ''<sub>''n'' +1                            </sub> (Ω<sup>(''e'' )                            </sup> ) representa el dominio deformado actual. Para detalles de linealización pueden consultarse textos estándar para deformaciones no lineales como Bonet y Wood [[#bib0110|[22]]] . El método de Newton-Raphson se usa como técnica implícita también en este caso conteniendo anidado el algoritmo de retorno [[#bib0115|[23]]]  (deformaciones finitas) para la resolución del sistema de ecuaciones de la forma débil del problema estructural:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{K}}_T\delta \boldsymbol{u}^{\left(k\right)}=</math><math>-\boldsymbol{R}^{\left(k-1\right)}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)

|}

​

donde '''K'''<sub>''T''</sub>  es obtenido en la ecuación [[#eq0190|(36)]]  a través de '''G'''  (operador espacial discreto). En el caso de tensión o deformación plana, '''G'''  toma la forma representada en la ecuación [[#eq0195|(37)]] :

​

<span id='eq0190'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{K}}_T\quad ={\bigwedge}_{e=1}^{nelem}\left\{{\int }_{{\varphi }_{n+1}({\Omega }^{\left(e\right)})}\boldsymbol{\mbox{G}}^T\boldsymbol{\mbox{a}}\boldsymbol{\mbox{G}}dv\right\}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)

|}

​

<span id='eq0195'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{G}}=\left[\begin{array}{ccccccc}

N_{1,1}^{\left(e\right)} & 0 & N_{2,1}^{\left(e\right)} & 0 & \ldots  & N_{n_{node},1}^{\left(e\right)} & 0\\

0 & N_{1,1}^{\left(e\right)} & 0 & N_{2,1}^{\left(e\right)} & \ldots  & 0 & N_{n_{node},1}^{\left(e\right)}\\

N_{1,2}^{\left(e\right)} & 0 & N_{2,2}^{\left(e\right)} & 0 & \ldots  & N_{n_{node},2}^{\left(e\right)} & 0\\

0 & N_{1,2}^{\left(e\right)} & 0 & N_{2,2}^{\left(e\right)} & \ldots  & 0 & N_{n_{node},2}^{\left(e\right)}\\

 &  &  &  &  &  & 

\end{array}\right]</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)

|}

​

​

[[Image:draft_Content_721884677-1-s2.0-S0213131513000217-fx2.jpg|center|367px|Image for unlabelled figure]]

​

El tensor de cuarto orden '''a'''  es el módulo tangente consistente, el cual se define en coordenadas cartesianas en la ecuación [[#eq0200|(38)]]  al final de la iteración (''k''  − 1):

​

<span id='eq0200'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>a_{ijkl}=\frac{1}{J}\frac{\partial {\tau }_{ij}}{\partial F_{km}}F_{lm}-</math><math>{\sigma }_{il}{\delta }_{jk}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)

|}

​

Se presentan en el cuadro II las principales modificaciones en el caso de grandes deformaciones.      

​

==8. Algoritmo explícito==

​

La forma débil que resulta de un sistema de ecuaciones discretas de momento, correspondiente a la malla de elementos finitos, es a su vez discretizada en tiempo mediante diferencias finitas centradas en el instante genérico ''t''<sub>''n''</sub> :

​

<span id='eq0205'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{M}}\quad \ddot{\boldsymbol{\mbox{u}}}(t_n)+</math><math>\boldsymbol{\mbox{C}}\quad \dot{\boldsymbol{\mbox{u}}}(t_n)+</math><math>\boldsymbol{\mbox{f}}^{int}(\boldsymbol{\mbox{u}}_n)=</math><math>\boldsymbol{\mbox{f}}^{ext}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)

|}

​

donde '''M'''  es la matriz (diagonal) de masa, '''C'''  la matriz de amortiguamiento, '''f'''<sup>''int''</sup> ('''u'''<sub>''n''</sub> ) el vector de fuerzas internas nodales, '''f'''<sup>''ext''</sup>  el vector de fuerzas externas aplicadas en los nodos, y <math display="inline">\ddot{\boldsymbol{\mbox{u}}}(t_n),\dot{\boldsymbol{\mbox{u}}}(t_n),\boldsymbol{\mbox{u}}_n</math>  son, respectivamente, aceleraciones, velocidades y desplazamientos en los nodos:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\dot{\boldsymbol{\mbox{u}}}(t_{n-1/2})=\frac{\boldsymbol{\mbox{u}}(t_n)-\boldsymbol{\mbox{u}}(t_{n-1})}{\Delta t_n}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)

|}

​

Las ecuaciones (desacopladas) que deben resolverse en cada grado de libertad son:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\dot{u}}_{i,n+1/2}=(\frac{m_{ii}}{\Delta t_n}+\frac{C_{ii}}{2})^{-1}\left(\left(\frac{m_{ii}}{\Delta t_n}-\right. \right. </math><math>\left. \left. \frac{C_{ii}}{2}){\dot{u}}_{i,n-1/2}+\right. \right. </math><math>\left. \left. f_i^{ext}(t_n)-f_i^{int}(t_n\right)\right)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

Las condiciones de contorno son fácilmente aplicadas, mientras que las condiciones iniciales provienen de la última iteración que convergió en el flujo implícito. Las condiciones iniciales correspondientes a la transferencia desde el flujo implícito al flujo explícito son descritas en la sección [[#sec0090|9]] .      

​

===8.1. Matriz de masa===

​

Se necesita una matriz de masa diagonal para obtener un sistema de ecuaciones que se pueda resolver explícitamente. Lo mismo se aplica a la matriz de amortiguamiento. La matriz de masa se obtiene como:

​

<span id='eq0220'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{M}}={\bigwedge}_{e=1}^{nelem}{\int }_{{\Omega }_{\left(e\right)}}{\rho }_0\boldsymbol{\mbox{N}}_{\left(e\right)}^T\boldsymbol{\mbox{N}}_{\left(e\right)}\quad d\Omega </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)

|}

​

La matriz obtenida mediante la ecuación [[#eq0220|(41)]]  no es diagonal y, por tanto, se necesita una aproximación para hacerla diagonal. En este caso se ha optado por la técnica que permite obtener la componente de la diagonal como la suma de las componentes de la misma fila (véase Belytschko et al. [[#bib0120|[24]]] ).      

​

===8.2. Disipación espuria===

​

La disipación espuria asociada a las oscilaciones producidas en la resolución explícita temporal de las ecuaciones de momento dinámicas puede afectar sensiblemente a la conservación de la energía. Dichas oscilaciones espurias producen un error numérico que puede reducirse mediante el uso de un amortiguamiento adecuado dependiendo del rango de frecuencias naturales de la estructura. Para reducir al mínimo la disipación espuria y obtener la solución del problema cuasi-estático se han considerado diversas técnicas de amortiguamiento, descritas en la siguiente sección. Se usa la técnica de relajación dinámica (Underwood [[#bib0125|[25]]] ) para la obtención de la respuesta estacionaria (correspondiente a la solución del problema cuasi-estático) mediante el uso de amortiguamiento en la solución del problema dinámico. Otras técnicas se basan en la introducción de una viscosidad artificial para la reducción de la susodicha disipación. Así, por ejemplo, el uso de viscosidad artificial ha sido satisfactoriamente empleado en la solución de problemas de impacto a baja velocidad (Nsiampa et al. [[#bib0130|[26]]] ).      

​

===8.3. Matriz de amortiguamiento===

​

La integración explícita de las ecuaciones de momento puede presentar extrañas oscilaciones que pueden, no obstante, ser controladas mediante el amortiguamiento. Una discusión sobre la estabilización de métodos numéricos mediante un amortiguamiento artificial puede verse en Park [[#bib0135|[27]]] . Además, en el caso de materiales no lineales es preciso usar un amortiguamiento que pueda variar con la rigidez, como el que señalan Belytschko et al. [[#bib0140|[28]]] . La matriz de amortiguamiento '''C'''  puede ser elegida diagonal para obtener una ventaja computacional y seguir teniendo un sistema de ecuaciones desacoplado que permita la solución explícita:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{C}}=\alpha \boldsymbol{\mbox{M}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)

|}

​

En concreto ''α'' , en el caso de matrices diagonales, puede ser modelado como:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\alpha }_i=2\xi {\omega }_i</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)

|}

​

donde ''ξ''  es la tasa de amortiguamiento y ''ω''<sub>''i''</sub>  es dado como:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\omega }_i^2=K_i/M_i</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)

|}

​

Existen otras opciones, como la propuesta por Munjiza et al. [[#bib0145|[29]]] , que considera el efecto de la rigidez:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{C}}=\boldsymbol{\mbox{M}}{\left(\boldsymbol{\mbox{M}}^{-1}\boldsymbol{\mbox{K}}\right)}^m</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)

|}

​

Para obtener la solución estática, usando la ecuación de momento dinámica [[#eq0205|(39)]] , se ha usado la relajación dinámica [[#bib0150|[30]]]  en la cual el estado estacionario después de un estado transitorio inicial corresponde a la solución estática. Así, si ''m''  = 1, se produce un amortiguamiento proporcional a la rigidez, lo cual amortigua las altas frecuencias. Sin embargo, el paso de tiempo crítico disminuye a medida que aumenta el amortiguamiento, lo que incrementa el coste computacional. Simulaciones con ''m''  = 0.5 produjeron en general un amortiguamiento de un amplio rango de frecuencias. ''m''  puede ser adaptado a los requerimientos del problema particular.      

​

===8.4. Paso de tiempo===

​

El paso de tiempo en el método explícito tiene que ser menor que un valor crítico para garantizar la estabilidad y, por tanto, la convergencia. Este valor se define en función de las frecuencias naturales ''ω''<sub>''i''</sub>  y de la tasa de amortiguamiento ''ξ''<sub>''i''</sub>  [véase ecuación [[#eq0245|(46)]] ], en cada nodo ''i''  de la malla:

​

<span id='eq0245'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta t\leq min\frac{2}{{\omega }_i}(-{\xi }_i+\sqrt{1+{\xi }_i^2})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 46)

|}

​

Esta inecuación se satisface si se elige la máxima frecuencia de la malla, que se puede calcular conociendo el máximo autovalor de la malla <math display="inline">{\omega }_{max}=\sqrt{{\lambda }_{max}^{\left(mesh\right)}}</math> . Además, el máximo autovalor puede ser delimitado por el máximo autovalor de los elementos de la malla <math display="inline">{\lambda }_{max}^{\left(mesh\right)}\leq {\lambda }_{max}^{\left(e\right)}</math>  (véase Irons y Ahmad [[#bib0155|[31]]]  o Irons y Treharne [[#bib0160|[32]]] ). El problema de autovalores añade más coste computacional. Para evitar elevar el coste, se ha llevado a cabo una estrategia alternativa. Para problemas elásticos se puede usar:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta t\leq min\frac{L_e}{c_e}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 47)

|}

​

donde,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>c_e=\sqrt{\frac{E^{\left(e\right)}}{{\rho }^{\left(e\right)}}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 48)

|}

​

donde ''L''<sub>''e''</sub>  es una longitud característica del elemento y ''c''<sub>''e''</sub>  es la velocidad de onda. Dicho paso de tiempo es denominado en la literatura ''número de Courant'' . Para el caso de deformación plana:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta t\leq L_e\sqrt{\frac{\rho (1+\nu )(1-2\nu )}{E(1-\nu )}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 49)

|}

​

y para tensión plana,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta t\leq L_e\sqrt{\frac{\rho (1-{\nu }^2)}{E}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 50)

|}

​

Estos pasos de tiempo son constantes y por tanto no recogen los nuevos cambios de configuración ni de material, por lo que es recomendable la adaptación del paso de tiempo crítico. Así, una posibilidad es:

​

<span id='eq0270'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta t\quad (t_{n+1})=\frac{2}{{max}_i\lbrace {\omega }_i(t_n)\rbrace }</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 51)

|}

​

Las frecuencias naturales se determinan resolviendo el problema homogéneo [ecuación [[#eq0275|(52)]] ]. Su solución analítica es de la forma <math display="inline">u(t)=\overset{\texttildelow ;}{u}e^{-j\omega t}</math><math display="inline">(j=\sqrt{-1})</math> , la cual, sustituida en la ecuación [[#eq0275|(52)]] , provee el problema de autovalores nuevamente en la ecuación [[#eq0285|(54)]] :

​

<span id='eq0275'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{Mü}}+\boldsymbol{\mbox{ku}}=0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 52)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>-{\omega }^2\boldsymbol{\mbox{M}}(cos(\omega t)-jsin(\omega t))+</math><math>\boldsymbol{\mbox{K}}(cos(\omega t)-jsin(\omega t))=</math><math>0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 53)

|}

​

<span id='eq0285'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\vert -{\omega }^2\boldsymbol{\mbox{M}}+\boldsymbol{\mbox{K}}\vert =</math><math>0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 54)

|}

​

donde ''ω''<sup>2</sup>  son los autovalores del sistema que proveen las frecuencias asociadas a cada grado de libertad. Para evitar el problema de autovalores hemos aproximado las rigideces como sigue:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>K_{ii}(t_n)\simeq \frac{f_i^{int}(u_n)-f_i^{int}(u_{n-1})}{{\dot{u}}_i^{n-1/2}\Delta t_n}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 55)

|}

​

Así, las frecuencias naturales pueden ser calculadas sin coste adicional como:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\omega }_i(t_n)=\sqrt{\frac{K_i(t_n)}{M_i}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 56)

|}

​

El paso de tiempo crítico se determina sustituyendo las frecuencia natural máxima en la ecuación [[#eq0270|(51)]] .      

​

<span id='sec0090'></span>

==9. Transición entre flujos implícito y explícito==

​

La transición se puede producir en cualquier instante dentro de un determinado incremento de carga. Así, el método de Newton-Raphson (con algoritmo de retorno anidado en el bucle de actualización de variables internas) podría diverger debido a que la norma del residual rebasa la tolerancia permitida en 2 iteraciones consecutivas y, como consecuencia, el método implícito pararía. En este caso, el flujo de procesamiento es desviado hacia la resolución explícita temporal de las ecuaciones de momento dinámicas. Para obtener la respuesta estacionaria correspondiente al problema cuasi-estático se ha usado la relajación dinámica [[#bib0125|[25]]] . El método explícito comienza con condiciones iniciales correspondientes al último incremento de carga que convergió mediante el método implícito. Así, las últimas iteraciones de Newton-Raphson que divergieron no son tenidas en cuenta en la inicialización, que comienza de esta forma desde un punto estable, es decir, los valores nodales de la última iteración que convergió en el flujo implícito (<math display="inline">\overset{\texttildelow ;}{\boldsymbol{\mbox{u}}}</math> ) son transferidos al flujo explícito como condiciones iniciales. Esto significa que los desplazamientos y las fuerzas internas en la última iteración del implícito no están en equilibrio. Las fuerzas internas nodales y los desplazamientos desde el implícito son usados como sigue:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{f}}^{int}(\overset{\texttildelow ;}{\boldsymbol{\mbox{u}}})\rightarrow \boldsymbol{\mbox{f}}^{int}(0)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 57)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\overset{\texttildelow ;}{\boldsymbol{\mbox{u}}}\rightarrow \boldsymbol{\mbox{u}}(0)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 58)

|}

​

Los desplazamientos y las fuerzas internas nodales correspondientes al susodicho punto estable son utilizados en la computación de las aceleraciones y velocidades nodales para el inicio del flujo explícito como sigue:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\ddot{u}}_i(0)=\frac{f_i^{ext}-f_i^{int}(0)}{M_{ii}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 59)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\dot{u}}_i(0)={\ddot{u}}_i(0)\Delta t(0)+{\dot{u}}_i^{-}(0)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 60)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\dot{u}}_i^{-}(0)=0.0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 61)

|}

​

donde Δ''t'' (0) es el paso de tiempo inicial. Después de la primera iteración, se lleva a cabo un paso de tiempo adaptativo.      

​

Una vez que la solución se ha alcanzado para el incremento de carga actual, el flujo es revertido de nuevo a la técnica implícita más rápida para el tipo de problemas considerado. En el caso de que la carga externa haya sido totalmente aplicada, se sale del bucle de carga y se finaliza la ejecución.

​

==10. Resultados numéricos==

​

Un programa de elementos finitos que integra la técnica propuesta ha sido codificado usando el lenguaje de programación FORTRAN. Para el post-proceso de resultados se creó una interface también en FORTRAN que volcara los resultados a ''GiD''[[#bib0165|[33]]]  para su visualización. Se presentan a continuación varios test numéricos en detalle.      

​

===10.1. Test 1: Problema de ''snap-trough''===

​

Para validar la técnica, se empieza por simular ''snap-through''  con un material hiperelástico. En este problema se presenta un arco esbelto usando un modelo de material neohookeano (módulo de Young ''E''  = 6.895 · 10<sup>4</sup>  ''MPa'' , coeficiente de Poisson ''ν''  = 0.34 y densidad ''ρ''  = 2.700 ''kg'' /''m''<sup>3</sup> ). La carga externa puntual es aplicada en el punto medio del mismo, como puede verse en la [[#fig0010|figura 2]] , hasta alcanzar un valor máximo ''F''  = 4.000 ''N'' . Algunos de los puntos de la curva de equilibrio se pueden obtener por supuesto sin problemas, como sucede en la respuesta lineal elástica inicial. Sin embargo, en los puntos límites Newton-Raphson divergió y la solución se obtuvo con el algoritmo combinado. El arco está totalmente anclado en los 2 soportes, sin posibilidad de rotación. Los parámetros geométricos de la sección son ''A''  = 25, 806 · 10<sup>−4</sup>  ''m''<sup>2</sup> , momento de inercia ''I''  = 5.54 · 10<sup>−7</sup>  ''m''<sup>4</sup> , espesor ''t''  = 0.0508 ''m'' , y radio de curvatura ''R''  = 5.08 ''m'' .

​

<span id='fig0010'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_721884677-1-s2.0-S0213131513000217-gr2.jpg|center|204px|Geometría del arco compuesto por 32 cuadriláteros con 4 puntos de integración.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 2.

​

Geometría del arco compuesto por 32 cuadriláteros con 4 puntos de integración.

​

</span>

|}

​

Este mismo problema ha sido estudiado por Calhoun y DaDeppo [[#bib0170|[34]]]  o Wen y Suhendro [[#bib0175|[35]]] , aunque solo proporcionaron el comportamiento hasta el punto límite. Un estudio más avanzado fue realizado por Pin y Trahair [[#bib0180|[36]]] , proporcionando puntos de la curva de equilibrio después de alcanzar el punto límite. En la simulación con la técnica propuesta se pudo constatar cómo Newton-Raphson diverge cuando la deflexión en el centro del arco es |''δ''<sub>''crit''</sub> | = 0.076 ''m  ''  correspondiente a una fuerza interna en dirección vertical <math display="inline">\boldsymbol{\mbox{f}}_{central\quad node,y}^{int}=</math><math>2781.917N</math> . En este punto, el flujo explícito es iniciado como indicado en la sección [[#sec0090|9]] . El resultado es mostrado en la [[#fig0015|figura 3]] .

​

<span id='fig0015'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_721884677-1-s2.0-S0213131513000217-gr3.jpg|center|374px|En esta figura se puede observar la evolución de las iteraciones para la ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 3.

​

En esta figura se puede observar la evolución de las iteraciones para la deflexión central del arco y la transición después de 3 iteraciones mediante Newton-Raphson. En la cuarta iteración, diverge y el flujo explícito es iniciado.

​

</span>

|}

​

El flujo explícito, iniciado después de 3 iteraciones del implícito, puede observarse en la [[#fig0020|figura 4]] . La solución para una carga de 4.000 ''N''  se alcanzó en un valor de la deflexión en el nodo central del arco igual a |''δ''<sub>''sol''</sub> | = 1.17 ''m'' , como puede verse en la  [[#fig0025|figura 5]] . Otras variables son presentadas en las figuras  [[#fig0035|Figura 7]] , [[#fig0040|Figura 8]]  and [[#fig0045|Figura 9]] . Además, varios tests se hicieron con cargas más pequeñas hasta obtener los puntos de la curva de equilibrio representada en la [[#fig0030|figura 6]] , donde se puede observar claramente el problema de ''snap-through'' . Los detalles de la ejecución del proceso pueden encontrarse en la  [[#tbl0005|tabla 1]] , la cual destaca el intercambio de técnicas implícita y explícita y los parámetros de convergencia. El tiempo de proceso del explícito fue de 21,32 segundos. El número de pasos del explícito fue de 42.644.