Source: https://www.scipedia.com/public/Curiel-Sosa_Owen_2013a
Timestamp: 2019-08-18 09:21:12+00:00

Document:
Switching implicit-explicit numerical technique for nonlinear unstable structures - Curiel-Sosa Owen 2013a - Scipedia
J.L. Curiel-Sosa, D.R.J. Owen, Conmutación temporal entre técnicas numéricas implícita y explícita para simulación de estructuras no lineales inestables, Rev. int. métodos numér. cálc. diseño ing., 29(2) (2013), p 92-103. URL https://www.scipedia.com/public/Curiel-Sosa_Owen_2013a
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.
Método de elementos finitos ; Inestabilidad ; Estructuras ; Implícito ; Explícito
Finite element method ; Instability ; Structures ; Implicit ; Explicit
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 [1] and [2] , los cuales conllevan con frecuencia problemas de estabilidad aparte de ser costosos. Puede consultarse en Kojić [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. [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 [5] lo aplicó con éxito a la resolución de problemas de reacción-difusión. Más recientemente, Robinson [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. [7] aplican una discretización temporal implícita-explícita usando elementos finitos para el problema de propagación de ondas. Comblen et al. [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. [9] o Farhat et al. [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 [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 [12] o de Laursen y Simo [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 [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 [15] , [16] , [17] and [18] . Matías y Oñate [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 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:
E = { σ ∥ Y ( σ , α ) < 0 } {\displaystyle E=\lbrace {\boldsymbol {\sigma }}\parallel Y({\boldsymbol {\sigma }},{\boldsymbol {\alpha }})\quad <0\rbrace }
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:
P = { σ ∥ Y ( σ , α ) = {\displaystyle P=\lbrace {\boldsymbol {\sigma }}\parallel Y({\boldsymbol {\sigma }},{\boldsymbol {\alpha }})=} 0 } {\displaystyle 0\rbrace }
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:
ϵ ˙ p = γ ˙ N ( σ , α ) = {\displaystyle {\boldsymbol {\dot {\epsilon }}}^{p}={\dot {\gamma }}\quad {\boldsymbol {\mbox{N}}}({\boldsymbol {\sigma }},{\boldsymbol {\alpha }})=} γ ˙ ∂ Ψ ∂ σ {\displaystyle {\dot {\gamma }}\quad {\frac {\partial \Psi }{\partial {\boldsymbol {\sigma }}}}}
donde N (σ , α ) es el vector de flujo. La ley de endurecimiento es:
α ˙ = γ ˙ H ( σ , α ) {\displaystyle {\boldsymbol {\dot {\alpha }}}={\dot {\gamma }}\quad {\boldsymbol {\mbox{H}}}({\boldsymbol {\sigma }},{\boldsymbol {\alpha }})}
donde H (σ , α ) define las variables de endurecimiento y es llamado módulo de endurecimiento generalizado [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 (3) , y de endurecimiento, ecuación (4) , se complementa con el criterio de carga/descarga, representado por las condiciones de Kuhn-Tucker:
Φ ≤ 0 γ ˙ ≥ 0 Φ γ ˙ = {\displaystyle \Phi \quad \leq \quad 0\quad {\dot {\gamma }}\geq \quad 0\quad \Phi \quad {\dot {\gamma }}=} 0 {\displaystyle 0}
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 J2 (s ) alcanza un valor crítico σy . 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, σ1 = − σ2 > 0, σ3 = 0 y J 2 ( s ) = {\textstyle {\sqrt {J_{2}({\boldsymbol {\mbox{s}}})}}=} τ m a x = σ 1 {\displaystyle {\tau }_{max}={\sigma }_{1}} la función toma la forma:
Φ ( σ ) = J 2 ( s ) − {\displaystyle \Phi ({\boldsymbol {\sigma }})={\sqrt {J_{2}({\boldsymbol {\mbox{s}}})}}-} τ y {\displaystyle {\tau }^{y}}
En el caso de tensiones unidireccionales, el criterio se postula:
Φ ( σ ) = 3 J 2 ( s ) − {\displaystyle \Phi ({\boldsymbol {\sigma }})={\sqrt {3\quad J_{2}({\boldsymbol {\mbox{s}}})}}-} σ y {\displaystyle {\sigma }^{y}}
donde σy es el valor del límite elástico obtenido en el ensayo de tracción. 3 J 2 ( s ) {\textstyle {\sqrt {3\quad J_{2}({\boldsymbol {\mbox{s}}})}}} es conocido como la tensión efectiva de Von-Mises o tensión equivalente. Véase que:
σ y = 3 τ y {\displaystyle {\sigma }^{y}={\sqrt {3}}\quad {\tau }^{y}}
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 [tn , tn +1 ], el incremento de la tasa de deformaciones es dado por la ecuación (8) . Las variables de estado consideradas son la deformación elástica ϵ n e {\textstyle {\boldsymbol {\epsilon }}_{n}^{e}} y la deformación plástica acumulada ϵ ¯ n p {\textstyle {\boldsymbol {\overline {\epsilon }}}_{n}^{p}} al principio del intervalo de tiempo [tn , tn +1 ]. La predicción del estado elástico de deformaciones viene dada por la ecuación (9) . La ecuación (10) recoge que no se incremente la deformación plástica:
Δ ϵ = ϵ n + 1 − {\displaystyle \Delta {\boldsymbol {\epsilon }}={\boldsymbol {\epsilon }}_{n+1}-} ϵ n {\displaystyle {\boldsymbol {\epsilon }}_{n}}
ϵ n + 1 e t r i a l = ϵ n e + {\displaystyle {\boldsymbol {\epsilon }}_{n+1}^{etrial}={\boldsymbol {\epsilon }}_{n}^{e}+} Δ ϵ {\displaystyle \Delta {\boldsymbol {\epsilon }}}
ϵ ¯ n + 1 p t r i a l = {\displaystyle {\boldsymbol {\overline {\epsilon }}}_{n+1}^{ptrial}=} ϵ ¯ n p {\displaystyle {\boldsymbol {\overline {\epsilon }}}_{n}^{p}}
El tensor de tensiones predecido asumiendo una evolución plástica viene dado por la ecuación (11) . La tensión correspondiente al límite elástico depende de la deformación plástica acumulada en el instante tn , ecuación (12) :
σ n + 1 t r i a l = D e : ϵ n + 1 e t r i a l {\displaystyle {\boldsymbol {\sigma }}_{n+1}^{trial}={\boldsymbol {\mbox{D}}}^{e}:{\boldsymbol {\epsilon }}_{n+1}^{etrial}}
σ y n + 1 t r i a l = σ y ( ϵ ¯ n p ) {\displaystyle {\sigma }_{yn+1}^{trial}={\sigma }_{y}({\boldsymbol {\overline {\epsilon }}}_{n}^{p})}
Si σ n + 1 t r i a l {\textstyle {\boldsymbol {\sigma }}_{n+1}^{trial}} yace en el interior del dominio elástico, ecuación (13) , la evolución es puramente elástica en el intervalo de tiempo [tn , tn +1 ] y la predicción corresponde a la solución para ese incremento:
Φ ( σ n + 1 t r i a l , σ y n ) ≤ 0 {\displaystyle \Phi ({\boldsymbol {\sigma }}_{n+1}^{trial},{\sigma }_{yn})\leq 0}
La actualización es simplemente llevada a cabo como sigue:
ϵ n + 1 e = ϵ n + 1 e t r i a l {\displaystyle {\boldsymbol {\epsilon }}_{n+1}^{e}={\boldsymbol {\epsilon }}_{n+1}^{etrial}}
ϵ ¯ n + 1 p = ϵ ¯ n + 1 p t r i a l = {\displaystyle {\boldsymbol {\overline {\epsilon }}}_{n+1}^{p}={\boldsymbol {\overline {\epsilon }}}_{n+1}^{ptrial}=} ϵ ¯ n p {\displaystyle {\boldsymbol {\overline {\epsilon }}}_{n}^{p}}
σ n + 1 = σ n + 1 t r i a l {\displaystyle {\boldsymbol {\sigma }}_{n+1}={\boldsymbol {\sigma }}_{n+1}^{trial}}
σ y n + 1 = σ y n + 1 t r i a l = σ y n {\displaystyle {\sigma }_{yn+1}={\sigma }_{yn+1}^{trial}={\sigma }_{yn}}
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 fig. 1 ) es necesario, como sigue:
ϵ n + 1 e = ϵ n + 1 e t r i a l − {\displaystyle {\boldsymbol {\epsilon }}_{n+1}^{e}={\boldsymbol {\epsilon }}_{n+1}^{etrial}-} Δ γ 2 3 s n + 1 ∥ s n + 1 ∥ {\displaystyle \Delta \gamma {\sqrt {\frac {2}{3}}}{\frac {{\boldsymbol {\mbox{s}}}_{n+1}}{\parallel {\boldsymbol {\mbox{s}}}_{n+1}\parallel }}}
ϵ ¯ n + 1 p = ϵ ¯ n p + {\displaystyle {\boldsymbol {\overline {\epsilon }}}_{n+1}^{p}={\boldsymbol {\overline {\epsilon }}}_{n}^{p}+} Δ γ {\displaystyle \Delta \gamma }
0 = 3 J 2 ( s n + 1 ) − σ y ( ϵ ¯ n + 1 p ) {\displaystyle 0={\sqrt {3\quad J_{2}({\boldsymbol {\mbox{s}}}_{n+1})}}-{\sigma }_{y}({\boldsymbol {\overline {\epsilon }}}_{n+1}^{p})}
Pseudo integración para caracterización del endurecimiento: predicción elástica – corrección plástica.
Este grupo de ecuaciones algebraicas tiene que ser resuelto para ϵ n + 1 e {\textstyle {\boldsymbol {\epsilon }}_{n+1}^{e}} , ϵ ¯ n + 1 p {\textstyle {\boldsymbol {\overline {\epsilon }}}_{n+1}^{p}} y Δγ . sn +1 es el tensor desviador en tn +1 ecuación (21) :
s n + 1 = 2 G d e v [ ϵ n + 1 e ] {\displaystyle {\boldsymbol {\mbox{s}}}_{n+1}=2Gdev[{\boldsymbol {\epsilon }}_{n+1}^{e}]}
Este sistema puede ser reducido a una sola ecuación, según Souza [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 [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 [20] . Así, si en un incremento de tiempo [tn , tn +1 ] se da un conjunto de variables internas αn en el instante tn , el tensor de deformaciones ε (tn +1 ) debe determinar las tensiones σ (tn +1 ) y las variables internas solo a través del algoritmo de integración. Véase Simo y Hughes [2] para detalles de este tipo de algoritmos, es decir:
σ ( t n + 1 ) = σ ˆ ( α n , ϵ n + 1 ) {\displaystyle {\boldsymbol {\sigma }}(t_{n+1})={\boldsymbol {\overset {\mbox{ˆ}}{\sigma }}}({\boldsymbol {\alpha }}_{n},{\boldsymbol {\epsilon }}_{n+1})}
α ( t n + 1 ) = α ˆ ( α n , ϵ n + 1 ) {\displaystyle {\boldsymbol {\alpha }}(t_{n+1})={\boldsymbol {\overset {\mbox{ˆ}}{\alpha }}}({\boldsymbol {\alpha }}_{n},{\boldsymbol {\epsilon }}_{n+1})}
Después de la discretización mediante elementos finitos, el problema se centra en encontrar los desplazamientos nodales un +1 en el instante tn +1 , de modo que se satisfaga el sistema de ecuaciones incremental no lineal (24) :
R ( u n + 1 ) = f i n t ( u n + 1 ) − {\displaystyle {\boldsymbol {R}}({\boldsymbol {u}}_{n+1})={\boldsymbol {\mbox{f}}}^{int}({\boldsymbol {u}}_{n+1})-} f n + 1 e x t = 0 {\displaystyle {\boldsymbol {\mbox{f}}}_{n+1}^{ext}=0}
donde las fuerzas nodales internas fint (un +1 ) y externas f n + 1 e x t {\textstyle {\boldsymbol {\mbox{f}}}_{n+1}^{ext}} se obtienen como sigue:
f i n t ( u n + 1 ) = {\displaystyle {\boldsymbol {\mbox{f}}}^{int}({\boldsymbol {u}}_{n+1})=} ⋀ e = 1 n e l e m { ∫ Ω ( e ) B T σ ˆ ( α n , ϵ ( u n + 1 ) ) d v } {\displaystyle {\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\}}
f n + 1 e x t = ⋀ e = 1 n e l e m { ∫ Ω ( e ) N T b n + 1 d v + {\displaystyle {\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.} ∫ ∂ Ω ( e ) N T q n + 1 d s } {\displaystyle \left.{\int }_{\partial {\Omega }^{\left(e\right)}}{\boldsymbol {\mbox{N}}}^{T}\quad {\boldsymbol {q}}_{n+1}ds\right\}}
donde N contiene las N (ξ , η ) funciones de formas bilineales, bn +1 son las fuerzas volumétricas, qn +1 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):
B = [ N 1 , 1 ( e ) 0 N 2 , 1 ( e ) 0 … N n n o d e , 1 ( e ) 0 0 N 1 , 2 ( e ) 0 N 2 , 2 ( e ) … 0 N n n o d e , 2 ( e ) N 1 , 2 ( e ) N 1 , 1 ( e ) N 2 , 2 ( e ) N 2 , 1 ( e ) … N n n o d e , 2 ( e ) N n n o d e , 1 ( e ) ] {\displaystyle {\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]}
Para continuar con el procedimiento numérico, la ecuación (24) necesita ser linealizada. Pueden encontrarse en Souza et al. [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 (24) ]. En el caso de materiales no lineales, en cada iteración hay que resolver la versión linealizada [ecuación (24) ] para el cálculo de desplazamientos nodales δ''u(k ) :
K T δ u ( k ) = {\displaystyle {\boldsymbol {\mbox{K}}}_{T}\delta {\boldsymbol {u}}^{\left(k\right)}=} − R ( k − 1 ) ( u n + 1 ) {\displaystyle -{\boldsymbol {R}}^{\left(k-1\right)}({\boldsymbol {u}}_{n+1})}
donde KT es la matriz tangente de rigidez dada por:
K T = ∂ R ∂ u n + 1 | u n + 1 ( k − 1 ) {\displaystyle {\boldsymbol {\mbox{K}}}_{T}={\frac {\partial {\boldsymbol {R}}}{\partial \quad {\boldsymbol {u}}_{n+1}}}\vert _{{\boldsymbol {u}}_{n+1}^{\left(k-1\right)}}}
la cual es obtenida mediante ensamblaje de las matrices tangentes de rigidez de cada elemento:
k T ( e ) = ∫ Ω ( e ) B T D ˆ B d v {\displaystyle {\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}
donde D ˆ {\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{D}}}}} es la matriz consistente tangente de rigidez [20] :
D ˆ = ∂ σ ˆ ∂ ϵ n + 1 | ϵ n + 1 ( k − 1 ) {\displaystyle {\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)}}}
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 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 σ n + 1 = {\textstyle {\boldsymbol {\sigma }}_{n+1}=} σ ˆ ( α n , F n + 1 ) {\displaystyle {\boldsymbol {\overset {\mbox{ˆ}}{\sigma }}}({\boldsymbol {\alpha }}_{n},{\boldsymbol {\mbox{F}}}_{n+1})} , donde la parte derecha corresponde a la tensión obtenida mediante el algoritmo de retorno. Fn +1 es el gradiente de deformaciones calculado al final del intervalo [tn , tn +1 ]. Ahora los vectores que contienen las fuerzas nodales (internas y externas) están basados en la configuración deformada:
f i n t ( u n + 1 ) = {\displaystyle {\boldsymbol {\mbox{f}}}^{int}({\boldsymbol {u}}_{n+1})=}
⋀ e = 1 n e l e m { ∫ φ n + 1 ( Ω ( e ) ) B T σ ( α n , F ( u n + 1 ) ) d v } {\displaystyle {\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\}}
f n + 1 e x t = ⋀ e = 1 n e l e m ∫ φ n + 1 ( Ω ( e ) ) N T b n + 1 d v + {\displaystyle {\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+}
+ ∫ ∂ φ n + 1 ( Ω ( e ) ) N T q n + 1 d s {\displaystyle +{\int }_{\partial {\varphi }_{n+1}({\Omega }^{\left(e\right)})}{\boldsymbol {\mbox{N}}}^{T}{\boldsymbol {q}}_{n+1}ds}
donde φn +1 (Ω(e ) ) representa el dominio deformado actual. Para detalles de linealización pueden consultarse textos estándar para deformaciones no lineales como Bonet y Wood [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 [23] (deformaciones finitas) para la resolución del sistema de ecuaciones de la forma débil del problema estructural:
K T δ u ( k ) = {\displaystyle {\boldsymbol {\mbox{K}}}_{T}\delta {\boldsymbol {u}}^{\left(k\right)}=} − R ( k − 1 ) {\displaystyle -{\boldsymbol {R}}^{\left(k-1\right)}}
donde KT es obtenido en la ecuación (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 (37) :
K T = ⋀ e = 1 n e l e m { ∫ φ n + 1 ( Ω ( e ) ) G T a G d v } {\displaystyle {\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\}}
G = [ N 1 , 1 ( e ) 0 N 2 , 1 ( e ) 0 … N n n o d e , 1 ( e ) 0 0 N 1 , 1 ( e ) 0 N 2 , 1 ( e ) … 0 N n n o d e , 1 ( e ) N 1 , 2 ( e ) 0 N 2 , 2 ( e ) 0 … N n n o d e , 2 ( e ) 0 0 N 1 , 2 ( e ) 0 N 2 , 2 ( e ) … 0 N n n o d e , 2 ( e ) ] {\displaystyle {\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]}
El tensor de cuarto orden a es el módulo tangente consistente, el cual se define en coordenadas cartesianas en la ecuación (38) al final de la iteración (k − 1):
a i j k l = 1 J ∂ τ i j ∂ F k m F l m − {\displaystyle a_{ijkl}={\frac {1}{J}}{\frac {\partial {\tau }_{ij}}{\partial F_{km}}}F_{lm}-} σ i l δ j k {\displaystyle {\sigma }_{il}{\delta }_{jk}}
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 tn :
M u ¨ ( t n ) + {\displaystyle {\boldsymbol {\mbox{M}}}\quad {\ddot {\boldsymbol {\mbox{u}}}}(t_{n})+} C u ˙ ( t n ) + {\displaystyle {\boldsymbol {\mbox{C}}}\quad {\dot {\boldsymbol {\mbox{u}}}}(t_{n})+} f i n t ( u n ) = {\displaystyle {\boldsymbol {\mbox{f}}}^{int}({\boldsymbol {\mbox{u}}}_{n})=} f e x t {\displaystyle {\boldsymbol {\mbox{f}}}^{ext}}
donde M es la matriz (diagonal) de masa, C la matriz de amortiguamiento, fint (un ) el vector de fuerzas internas nodales, fext el vector de fuerzas externas aplicadas en los nodos, y u ¨ ( t n ) , u ˙ ( t n ) , u n {\textstyle {\ddot {\boldsymbol {\mbox{u}}}}(t_{n}),{\dot {\boldsymbol {\mbox{u}}}}(t_{n}),{\boldsymbol {\mbox{u}}}_{n}} son, respectivamente, aceleraciones, velocidades y desplazamientos en los nodos:
u ˙ ( t n − 1 / 2 ) = u ( t n ) − u ( t n − 1 ) Δ t n {\displaystyle {\dot {\boldsymbol {\mbox{u}}}}(t_{n-1/2})={\frac {{\boldsymbol {\mbox{u}}}(t_{n})-{\boldsymbol {\mbox{u}}}(t_{n-1})}{\Delta t_{n}}}}
Las ecuaciones (desacopladas) que deben resolverse en cada grado de libertad son:
u ˙ i , n + 1 / 2 = ( m i i Δ t n + C i i 2 ) − 1 ( ( m i i Δ t n − {\displaystyle {\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.} C i i 2 ) u ˙ i , n − 1 / 2 + {\displaystyle \left.\left.{\frac {C_{ii}}{2}}){\dot {u}}_{i,n-1/2}+\right.\right.} f i e x t ( t n ) − f i i n t ( t n ) ) {\displaystyle \left.\left.f_{i}^{ext}(t_{n})-f_{i}^{int}(t_{n}\right)\right)}
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 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:
M = ⋀ e = 1 n e l e m ∫ Ω ( e ) ρ 0 N ( e ) T N ( e ) d Ω {\displaystyle {\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 }
La matriz obtenida mediante la ecuación (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. [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 [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. [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 [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. [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:
C = α M {\displaystyle {\boldsymbol {\mbox{C}}}=\alpha {\boldsymbol {\mbox{M}}}}
En concreto α , en el caso de matrices diagonales, puede ser modelado como:
α i = 2 ξ ω i {\displaystyle {\alpha }_{i}=2\xi {\omega }_{i}}
donde ξ es la tasa de amortiguamiento y ωi es dado como:
ω i 2 = K i / M i {\displaystyle {\omega }_{i}^{2}=K_{i}/M_{i}}
Existen otras opciones, como la propuesta por Munjiza et al. [29] , que considera el efecto de la rigidez:
C = M ( M − 1 K ) m {\displaystyle {\boldsymbol {\mbox{C}}}={\boldsymbol {\mbox{M}}}{\left({\boldsymbol {\mbox{M}}}^{-1}{\boldsymbol {\mbox{K}}}\right)}^{m}}
Para obtener la solución estática, usando la ecuación de momento dinámica (39) , se ha usado la relajación dinámica [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 ωi y de la tasa de amortiguamiento ξi [véase ecuación (46) ], en cada nodo i de la malla:
Δ t ≤ m i n 2 ω i ( − ξ i + 1 + ξ i 2 ) {\displaystyle \Delta t\leq min{\frac {2}{{\omega }_{i}}}(-{\xi }_{i}+{\sqrt {1+{\xi }_{i}^{2}}})}
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 ω m a x = λ m a x ( m e s h ) {\textstyle {\omega }_{max}={\sqrt {{\lambda }_{max}^{\left(mesh\right)}}}} . Además, el máximo autovalor puede ser delimitado por el máximo autovalor de los elementos de la malla λ m a x ( m e s h ) ≤ λ m a x ( e ) {\textstyle {\lambda }_{max}^{\left(mesh\right)}\leq {\lambda }_{max}^{\left(e\right)}} (véase Irons y Ahmad [31] o Irons y Treharne [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:
Δ t ≤ m i n L e c e {\displaystyle \Delta t\leq min{\frac {L_{e}}{c_{e}}}}
c e = E ( e ) ρ ( e ) {\displaystyle c_{e}={\sqrt {\frac {E^{\left(e\right)}}{{\rho }^{\left(e\right)}}}}}
donde Le es una longitud característica del elemento y ce 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:
Δ t ≤ L e ρ ( 1 + ν ) ( 1 − 2 ν ) E ( 1 − ν ) {\displaystyle \Delta t\leq L_{e}{\sqrt {\frac {\rho (1+\nu )(1-2\nu )}{E(1-\nu )}}}}
y para tensión plana,
Δ t ≤ L e ρ ( 1 − ν 2 ) E {\displaystyle \Delta t\leq L_{e}{\sqrt {\frac {\rho (1-{\nu }^{2})}{E}}}}
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:
Δ t ( t n + 1 ) = 2 m a x i { ω i ( t n ) } {\displaystyle \Delta t\quad (t_{n+1})={\frac {2}{{max}_{i}\lbrace {\omega }_{i}(t_{n})\rbrace }}}
Las frecuencias naturales se determinan resolviendo el problema homogéneo [ecuación (52) ]. Su solución analítica es de la forma u ( t ) = u ~ e − j ω t {\textstyle u(t)={\tilde {u}}e^{-j\omega t}} ( j = − 1 ) {\textstyle (j={\sqrt {-1}})} , la cual, sustituida en la ecuación (52) , provee el problema de autovalores nuevamente en la ecuación (54) :
Mü + ku = 0 {\displaystyle {\boldsymbol {\mbox{Mü}}}+{\boldsymbol {\mbox{ku}}}=0}
− ω 2 M ( c o s ( ω t ) − j s i n ( ω t ) ) + {\displaystyle -{\omega }^{2}{\boldsymbol {\mbox{M}}}(cos(\omega t)-jsin(\omega t))+} K ( c o s ( ω t ) − j s i n ( ω t ) ) = {\displaystyle {\boldsymbol {\mbox{K}}}(cos(\omega t)-jsin(\omega t))=} 0 {\displaystyle 0}
| − ω 2 M + K | = {\displaystyle \vert -{\omega }^{2}{\boldsymbol {\mbox{M}}}+{\boldsymbol {\mbox{K}}}\vert =} 0 {\displaystyle 0}
donde ω2 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:
K i i ( t n ) ≃ f i i n t ( u n ) − f i i n t ( u n − 1 ) u ˙ i n − 1 / 2 Δ t n {\displaystyle 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}}}}
Así, las frecuencias naturales pueden ser calculadas sin coste adicional como:
ω i ( t n ) = K i ( t n ) M i {\displaystyle {\omega }_{i}(t_{n})={\sqrt {\frac {K_{i}(t_{n})}{M_{i}}}}}
El paso de tiempo crítico se determina sustituyendo las frecuencia natural máxima en la ecuación (51) .
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 [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 ( u ~ {\textstyle {\tilde {\boldsymbol {\mbox{u}}}}} ) 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:
f i n t ( u ~ ) → f i n t ( 0 ) {\displaystyle {\boldsymbol {\mbox{f}}}^{int}({\tilde {\boldsymbol {\mbox{u}}}})\rightarrow {\boldsymbol {\mbox{f}}}^{int}(0)}
u ~ → u ( 0 ) {\displaystyle {\tilde {\boldsymbol {\mbox{u}}}}\rightarrow {\boldsymbol {\mbox{u}}}(0)}
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:
u ¨ i ( 0 ) = f i e x t − f i i n t ( 0 ) M i i {\displaystyle {\ddot {u}}_{i}(0)={\frac {f_{i}^{ext}-f_{i}^{int}(0)}{M_{ii}}}}
u ˙ i ( 0 ) = u ¨ i ( 0 ) Δ t ( 0 ) + u ˙ i − ( 0 ) {\displaystyle {\dot {u}}_{i}(0)={\ddot {u}}_{i}(0)\Delta t(0)+{\dot {u}}_{i}^{-}(0)}
u ˙ i − ( 0 ) = 0.0 {\displaystyle {\dot {u}}_{i}^{-}(0)=0.0}
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[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 · 104 MPa , coeficiente de Poisson ν = 0.34 y densidad ρ = 2.700 kg /m3 ). La carga externa puntual es aplicada en el punto medio del mismo, como puede verse en la 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−4 m2 , momento de inercia I = 5.54 · 10−7 m4 , espesor t = 0.0508 m , y radio de curvatura R = 5.08 m .
Geometría del arco compuesto por 32 cuadriláteros con 4 puntos de integración.
Este mismo problema ha sido estudiado por Calhoun y DaDeppo [34] o Wen y Suhendro [35] , aunque solo proporcionaron el comportamiento hasta el punto límite. Un estudio más avanzado fue realizado por Pin y Trahair [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 |δcrit | = 0.076 m correspondiente a una fuerza interna en dirección vertical f c e n t r a l n o d e , y i n t = {\textstyle {\boldsymbol {\mbox{f}}}_{central\quad node,y}^{int}=} 2781.917 N {\displaystyle 2781.917N} . En este punto, el flujo explícito es iniciado como indicado en la sección 9 . El resultado es mostrado en la 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.
El flujo explícito, iniciado después de 3 iteraciones del implícito, puede observarse en la 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 |δsol | = 1.17 m , como puede verse en la figura 5 . Otras variables son presentadas en las figuras Figura 7 , Figura 8 and 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 figura 6 , donde se puede observar claramente el problema de snap-through . Los detalles de la ejecución del proceso pueden encontrarse en la 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.
Valor absoluto de la deflexión en el centro. El flujo explícito empieza cuando se han ejecutado 3 iteraciones del método implícito (véase fig 3 ).
Evolución de las deformaciones del arco con valores de la fuerza interna y deflexión en el nodo central entre paréntesis.
Resultado del campo de desplazamientos horizontal ux (m ) para carga de 4.000 N .
Resultado para el campo de velocidades horizontales v x ( m / s ) {\textstyle v_{x}\quad (m/s)} para carga de 4.000 N . Se puede observar que los valores son prácticamente nulos, como corresponde a la solución estacionaria del método explícito, lo cual provee la solución al problema estático.
Resultado para el campo de velocidades verticales para carga de 4.000 N .
Curva de equilibrio obtenida usando la técnica propuesta.
Tabla 1. Norma relativa del residual(%) (error(%) en la tabla), residual máximo, fuerza interna nodal en dirección vertical en el punto medio f y i n t {\textstyle {\boldsymbol {\mbox{f}}}_{y}^{int}} (N ) y deflección uy (m ). En la segunda columna el número de iteración (i ) es mostrado si corresponde al flujo implícito o el paso de tiempo es enseñado si el flujo es explícito (tn )
i /tn
Max. Resid.
f y i n t {\textstyle {\boldsymbol {\mbox{f}}}_{y}^{int}}
IMP 1 3.118,820 197.857,00 - - -
IMP 2 59,346 3.767,75 -7.220,6 -0,061 divergiendo
IMP 3 3.116,590 202.416,00 -4.015,2 -0,074 conmutación
EXP 0,0001 73,404 1.616,89 -3.477,9 -0,074 oscilando
EXP 0,001 64,204 785,12 -4.389,4 -0,075 oscilando
EXP 0,1 2,386 336,19 -3.971,2 -0,466 oscilando
EXP 0,5 0,092 29,74 -3.998,5 -1,193 convergiendo
EXP 0,641 0,033 9,98 -3.999,3 -1,17 convergiendo
EXP 1,0 0,007 2,016 -4.001,27 -1,17 convergiendo
EXP 1,5 10−3 0,293 -4.000,16 -1,17 convergiendo
EXP 2,0 4, 4 · 10−5 0,017 -4.000,01 -1,17 solución
10.2. Test 2: Arco elastoplástico
En este caso se aplica un material no lineal en concreto: el modelo de Von-Mises asociativo isotrópico presentado en la sección 4 a la geometría de arco circular enseñada en la figura 10 . No se permite rotación ni desplazamientos en los soportes. El arco está formado por 20 elementos cuadriláteros con 4 puntos de integración. Se considera sujeto a tensión plana (espesor 10 cm). Los parámetros del material son dados en la tabla 2 .
Dimensiones del arco elastoplástico.
Tabla 2. Parámetros materiales del arco elastoplástico.
Densidad 7, 860 · 10 fracción kg/cm3
Módulo de Young 210 ⋅ 10 5 N c m 2 {\textstyle 210\cdot {10}^{5}\quad {\frac {N}{{cm}^{2}}}}
Límite elástico 24500 N c m 2 {\textstyle 24500\quad {\frac {N}{{cm}^{2}}}}
Cuando el arco es cargado puntualmente en su punto central (correspondiente al nodo más alto) con una carga de 6.6 N el método implícito empieza a diverger en la cuarta iteración. En ese momento, el flujo explícito es automáticamente iniciado (véase tabla 3 para mayor detalle). Sobre el campo de desplazamientos representado sobre la configuración original, véase la figura 12 .
Tabla 3. Para el arco elastoplástico: norma relativa del residual(%) (error(%) en la tabla), residual máximo, fuerza interna nodal en dirección vertical en el punto medio f y i n t {\textstyle {\boldsymbol {\mbox{f}}}_{y}^{int}} (N ) y deflexión uy (m ). En la segunda columna el número de iteración (i ) es mostrado si corresponde al flujo implícito o el paso de tiempo es enseñado si el flujo es explícito (tn )
IMP 1 47,48 237,70 - - inicial
IMP 2 25,72 163,08 -7430,40 -0,605 -
IMP 3 1,479 7336,16 -29695,92 -1,078 divergiendo
IMP 4 862,37 0, 44 · 10+7 -62912,7 -1,119 conmutación
EXP 0,3943256 972567 5493540 -5542582 -1,625 oscilando
EXP 3,33 0,802 2633,40 -63371,6 -1,522 -
EXP 6,28 0,45 1885,83 -64114,2 -1,703 -
EXP 99,99 0,000188 0,0096 -65998,2 -1,043 convergiendo
EXP 150 0,000008 0,0406 -65999,97 -1,044 convergiendo
EXP 170,56 1, 0 · 10−6 0,008 -65999,99 -1,044 solución
Campo de desplazamientos verticales en el arco elastoplástico. Fuerza en el nodo central igual a 6.6 · 104 N .
La gráfica logarítmica correspondiente al error relativo, fig. 11 enseña claramente el momento del cambio de flujo para la carga de 6.6 · 104 N . Como puede observarse, se detecta el aumento del error y se lleva a cabo la conmutación del flujo (desde implícito a explícito) reduciendo el error y convergiendo a la solución. Para cargas menores no hubo problemas de convergencia con el implícito. El tiempo de procesamiento del explícito fue de 46,2 segundos, correspondiente a un número de pasos de 92.476.
Error relativo de la norma del residual (arco elastoplástico).
10.3. Test 3: Grandes deformaciones en la membrana de Cook
La membrana de Cook ha sido usada frecuentemente para el seguimiento de la convergencia debido a la incompresibilidad cuando se carga en cortadura (véanse por ejemplo los trabajos de Simo y Armero [37] o Andrade et al. [38] ). En principio no hay problemas en obtener la solución con Newton-Raphson siempre y cuando se permita la división de la carga externa en diversos incrementos. Sin embargo, para ilustrar las capacidades de la técnica propuesta se descarta esa posibilidad, por lo que la carga es aplicada totalmente desde el principio y, consecuentemente, el método implícito diverge y es cambiado automáticamente al flujo explícito. La comparación de ambas soluciones debería proveer el mismo resultado con una lógica discrepancia. Las dimensiones de la membrana de Cook se pueden ver en la figura 13 . El modelo no lineal de Ogden [39] se usa para modelar la membrana de Cook.
Membrana de Cook (dimensiones en mm). Geometría y malla de elementos finitos.
La membrana está totalmente restringida en desplazamientos en su parte izquierda, y en su parte derecha se aplica una carga distribuida de cortadura de valor 100 N (véase fig. 13 ). El módulo de cortadura es μ = 80.1938 y κ = 40.0942 × 104 . Se llevaron a cabo 2 simulaciones. La primera, con división de la carga externa de cortadura para evitar divergencia de Newton-Raphson, con la que se obtiene una solución que se puede comparar. Y la segunda, como se ha mencionado anteriormente, no permite la división de la carga externa, y se usa la técnica combinada. Así, se dieron 4 iteraciones del implícito antes de desviar la ejecución del programa hacia el flujo explícito alcanzando la solución con un error del 0.2%. La solución puede ser observada en comparación con la obtenida mediante el implícito y la división de carga en las figuras Figura 14 , Figura 15 , Figura 16 , Figura 17 , Figura 18 , Figura 19 , Figura 20 , Figura 21 and Figura 22 .
Desplazamiento ux (mm) obtenido con la técnica propuesta.
Desplazamiento vertical uy (mm) mediante Newton-Raphson y división de carga externa.
Desplazamiento vertical uy (mm) mediante la técnica propuesta.
Tensión de cortadura σxy (N /mm2 ) mediante Newton-Raphson y división de carga externa.
Tensión de cortadura σxy (N /mm2 ) mediante la técnica propuesta.
Tensión σxx (N /mm2 ) mediante Newton-Raphson y división de carga externa.
Tensión σxx (N /mm2 ) mediante la técnica propuesta.
Tensión σyy (N /mm2 ) mediante Newton-Raphson y división de carga externa.
Tensión σyy (N /mm2 ) mediante la técnica propuesta.
El tiempo de proceso del explícito fue de 637 segundos frente a 42,4 segundos de la estrategía implícita. Las soluciones obtenidas son corroboradas en la literatura (véase, por ejemplo, Andrade et al. [38] ). Simo y Armero [37] usaron control de la carga con incrementos ΔF = 0.5 desde 0 a 100N y obtuvieron un desplazamiento para el nodo de referencia (esquina de arriba a la derecha) de 6.9 mm, coincidiendo con el obtenido aquí.
En este artículo se ha presentado una técnica combinada implícita-explícita para la solución mediante el MEF de problemas no lineales materiales y/o geométricos sujetos a deformaciones finitas. Concretamente, se estudiaron problemas de puntos críticos. La técnica combina la rapidez del método implícito usado con la convergencia (condicional) del explícito, con lo que se consigue una ejecución óptima para los problemas de interés abordados. Se han presentado asimismo detalles del problema incremental y de la transición entre flujos implícito y explícito. Los test numéricos sobre inestabilidad estructural y grandes deformaciones demostraron la capacidad de la técnica para abordar este tipo de problemas. A pesar de que no se ha estudiado todavía, esta técnica podría ser conveniente para problemas de contacto no suave entre esquina y superficie cóncava que hacen diverger al método implícito pero para los que tampoco es deseable ejecutar explícitamente en su totalidad debido al coste computacional.
El primer autor está agradecido por el soporte recibido por Rockfield Software Ltd.
[1] M. Ortiz, J.C. Simo; Analysis of a new class of integration algorithms for elastoplastic constitutive relations; Int. J. Num. Meth. Engng., 23 (3) (1986), pp. 353–366
[2] J.C. Simo, T.J.R. Hughes; Computational Inelasticity; Springer-Verlag, New York (1998)
[3] M. Kojić; Stress integration procedures for inelastic material models within the finite element method; Appl. Mechs. Reviews, 55 (4) (2002), pp. 389–414
[4] U.M. Ascher, S.J. Ruuth, B.T. Wetton; Implicit-explicit methods for time-dependent partial differential equations; SIAM J. Num. Anal., 32 (1995), pp. 797–823
[5] S.J. Ruuth; Implicit-Explicit Methods for Reaction-Diffusion Problems in Pattern Formation; J. Math. Biol., 34 (1995), pp. 148–176
[6] M. Robinson; IMEX method convergence for a parabollic equation; J. Diff. Equations, 241 (2) (2007), pp. 225–236
[7] A.V. Idesman, M. Scmidt, J.R. Foley; Accurate 3-D finite element simulation of elastic wave propagation with the combination of explicit and implicit time-integration methods; Wave Motion. (2011) doi: 10.1016/j.wavemoti.2011.04.017
[8] R. Comblen, S. Blaise, V. Legat, J-F. Remacle, E. Deleersnijder, J. Lambrechts; A discontinuous finite element baroclinic marine model on unstructured prismatic meshes; Ocean Dynamics, 60 (2010), pp. 1395–1414
[9] T.J.R. Hughes, K.S. Pister, R.L. Taylor; Implicit-explicit finite elements in nonlinear transient analysis; Comp. Meths. App. Mechs. Engrg., 17-18 (Part I) (1979), pp. 159–182
[10] C. Farhat, M. Lessoinne, N. Maman; Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution; Int. J. Numer. Meth. Fluids, 21 (10) (1995), pp. 807–835
[11] J.L. Curiel Sosa, E.A. de Souza Neto, D.R.J. Owen; A combined implicit-explicit algorithm in time for non-linear finite element analysis; Commun. Numer. Methods Eng., 22 (2006), pp. 63–75
[12] E. Oñate, F.G. Flores; Advances in the formulation of the rotation-free basic shell triangle; Comput. Methods Appl. Mech. Engrg., 194 (2005), pp. 2406–2443
[13] T.A. Laursen, J.C. Simo; A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems; Int. J. Num. Meth. Engng., 36 (20) (1993), pp. 3451–3485
[14] C.A. Felippa; Transverse critical points by penalty springs; M. Nijhoff (Ed.), Proceedings of NUMETA, Swansea, Dordrecht, Holland (1987)
[15] Y.T. Feng, D. Peric, D.R.J. Owen; Determination of travel directions in path-following methods; Int. J. Math. Comput. Modell., 21 (7) (1995), pp. 43–59
[16] Y.T. Feng, D. Peric, D.R.J. Owen; On the sign of the determinant of the structural stiffness matrix for determination of loading increment in arc-length algorithms; Commun. Numer. Meth. Engng., 13 (1997), pp. 47–49
[17] Y.T. Feng, D. Peric, D.R.J. Owen; A new criterion for determination of initial loading parameter in arc-length methods; Comput. Struct., 58 (3) (1995), pp. 479–485
[18] E.A. de Souza Neto, Y.T. Feng; On the determination of the path direction for arc-length methods in the presence of bifurcations and ‘snap-backs’; Comput. Methods Appl. Mech. Engrg., 179 (1999), pp. 81–89
[19] W. Matias, E. Oñate; Análisis de la inestabilidad de estructuras por el método de desplazamiento crítico; Rev. Int. Mét. Num. Cálc. Dis. Ing., 13 (2) (1997), pp. 241–263
[20] E.A. de Souza Neto, D. Peric, DR.J. Owen; Computational Methods for Plasticity: Theory and Applications; Wiley, Chichester, Reino Unido (2008)
[21] E.A. de Souza Neto; A fast, one-equation integration algorithm for the Lemaitre ductile damage model; Commun. Number. methods Eng., 18 (2002), pp. 541–554
[22] J. Bonet, R.D. Wood; Nonlinear Continuum Mechanics for Finite Element Analysis; Cambridge University Press, Cambridge, Reino Unido (1997)
[23] J.C. Simo; Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory; Comp. Meths. App. Mechs. Engrg., 99 (1) (1992), pp. 61–112
[24] T. Belytschko, W.K. Liu, B. Moran; Nonlinear Finite Elements for Continua and Structures; John Wiley and Sons, Chichester, Reino Unido (2001)
[25] P. Underwood; Dynamic relaxation; Computer Methods for Transient Analysis (1983), pp. 245–285
[26] N. Nsiampa, J. Ponthot, L. Noels; Comparative study of numerical explicit schemes for impact problems; Int. J. Impact Eng., 35 (12) (2008), pp. 1688–1694
[27] K.C. Park; Practical aspects of numerical time integration; Comput. Struct., 7 (1977), pp. 343–353
[28] T. Belytschko, R.L. Chiappetta, H.D. Bartel; Efficient large scale non-linear transient analysis by finite elements; Int. J. Numer. Meth. Engng., 10 (1976), pp. 579–596
[29] A. Munjiza, D.R.J. Owen, A.J.L. Crook; A M (M−1K )m proportional damping in explicit integration of dynamic structural systems ; Int. J. Num. Meth. Engng., 41 (7) (1998), pp. 1277–1296
[30] T. Belytschko, T.J.R. Hughes; Computational Methods for Transient Analysis; Elsevier Science B.V, The Netherlands (2001)
[31] B.M. Irons, S. Ahmad; Techniques of Finite Elements; Ellis Horwood Ltd. (1980)
[32] B.M. Irons, C. Treharne; A bound theorem for eigenvalues and its practical application; Proc. 3rd Conference of Matrix Methods in Structural Mechanics AFFDL-TR-71-160 Wright-Patterson, Ohio (1972), pp. 245–254
[33] GiD. The personal pre and post processor program. Disponible en: http://www.gid.cimne.upc.es [consultado 1 Feb 2009].
[34] P.R. Calhoun, D.A. DaDeppo; Nonlinear finite element analysis of clamped arches; J. Struct. Engng. ASCE, 109 (3) (1983), pp. 599–612
[35] R.K. Wen, B. Suhendro; Nonlinear curved-beam element for arch structures; J. Struct. Engng. ASCE, 117 (11) (1991), pp. 3496–3515
[36] Y-L. Pin, N.S. Trahair; Nonlinear buckling and postbuckling of elastic arches; Eng. Struct., 20 (7) (1998), pp. 571–579
[37] J.C. Simo, F. Armero; Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes; Int. J. Num. Meth. Engng., 33 (7) (1992), pp. 1413–1449
[38] F.M. Andrade Pires, E.A. de Souza Neto, J.L. de la Cuesta Padilla; An assessment of the average nodal volume formulation for the analysis of nearly incompressible solids under finite strain; Commun. Numer. Meth. Engng., 20 (2004), pp. 569–583
[39] R.W. Ogden; Large deformation isotropic elasticity - on the correlation of theory and experiment for incompressible rubberlike solids; Proc. R. Soc. Lond. A., 326 (1972), pp. 565–584
Accepted on 10/12/11
Submitted on 21/07/11
DOI: 10.1016/j.rimni.2013.04.006
Finite element method • Instability • Structures • Implicit • Explicit

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución