Conmutación temporal entre técnicas numéricas implícita y explícita para simulación de estructuras no lineales inestables

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.

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 [1,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]:

• 1.

  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.

• 2.

  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–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:

• 1.

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

• 2.

  Evolución elástica.

• 3.

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

• 4.

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

• 5.

  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:

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:

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:

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

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:

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₂(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, σ₁ = − σ₂ > 0, σ₃ = 0 y J2(s)=τmax=σ1 la función toma la forma:

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

donde σ^y es el valor del límite elástico obtenido en el ensayo de tracción. 3J2(s) es conocido como la tensión efectiva de Von-Mises o tensión equivalente. Véase que:

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_n, t_n+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 εne y la deformación plástica acumulada ε¯np al principio del intervalo de tiempo [t_n, t_n+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:

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 t_n, ecuación (12):

Si σn+1trial yace en el interior del dominio elástico, ecuación (13), la evolución es puramente elástica en el intervalo de tiempo [t_n, t_n+1] y la predicción corresponde a la solución para ese incremento:

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

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:

Figura 1.

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+1e, ε¯n+1p y Δγ. s_n+1 es el tensor desviador en t_n+1 ecuación (21):

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 [t_n, t_n+1] se da un conjunto de variables internas α_n en el instante t_n, el tensor de deformaciones ε(t_n+1) debe determinar las tensiones σ(t_n+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:

Después de la discretización mediante elementos finitos, el problema se centra en encontrar los desplazamientos nodales u_n+1 en el instante t_n+1, de modo que se satisfaga el sistema de ecuaciones incremental no lineal (24):

donde las fuerzas nodales internas f^int(u_n+1) y externas fn+1ext se obtienen como sigue:

donde N contiene las N(ξ, η) funciones de formas bilineales, b_n+1 son las fuerzas volumétricas, q_n+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):

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):

donde K_T es la matriz tangente de rigidez dada por:

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

donde Dˆ es la matriz consistente tangente de rigidez [20]:

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=σˆ(αn,Fn+1), donde la parte derecha corresponde a la tensión obtenida mediante el algoritmo de retorno. F_n+1 es el gradiente de deformaciones calculado al final del intervalo [t_n, t_n+1]. Ahora los vectores que contienen las fuerzas nodales (internas y externas) están basados en la configuración deformada:

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:

donde K_T 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):

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):

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_n:

donde M es la matriz (diagonal) de masa, C la matriz de amortiguamiento, f^int(u_n) el vector de fuerzas internas nodales, f^ext el vector de fuerzas externas aplicadas en los nodos, y u¨(tn),u˙(tn),un son, respectivamente, aceleraciones, velocidades y desplazamientos en los nodos:

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

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:

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:

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

donde ξ es la tasa de amortiguamiento y ω_i es dado como:

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

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:

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 ωmax=λmax(mesh). Además, el máximo autovalor puede ser delimitado por el máximo autovalor de los elementos de la malla λmax(mesh)≤λmax(e) (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:

donde,

donde L_e es una longitud característica del elemento y c_e 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:

y para tensión plana,

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:

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(j=−1), la cual, sustituida en la ecuación (52), provee el problema de autovalores nuevamente en la ecuación (54):

donde ω² 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:

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

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˜;) 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:

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:

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 · 10⁴ MPa, coeficiente de Poisson ν = 0.34 y densidad ρ = 2.700 kg/m³). 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⁻⁴ m², momento de inercia I = 5.54 · 10⁻⁷ m⁴, espesor t = 0.0508 m, y radio de curvatura R = 5.08 m.

Figura 2.

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 fcentralnode,yint=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.

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 7–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.

Figura 4.

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

Figura 5.

Evolución de las deformaciones del arco con valores de la fuerza interna y deflexión en el nodo central entre paréntesis.

Figura 7.

Resultado del campo de desplazamientos horizontal u_x (m) para carga de 4.000 N.

Figura 8.

Resultado para el campo de velocidades horizontales vx(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.

Figura 9.

Resultado para el campo de velocidades verticales para carga de 4.000 N.

Figura 6.

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 fyint (N) y deflección u_y (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 (t_n)

Flujo	i/tn	Error(%)	Max. Resid.	fyint	uy	Estado
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.

Figura 10.

Dimensiones del arco elastoplástico.

Tabla 2.

Parámetros materiales del arco elastoplástico.

Propiedades materiales
Densidad	7, 860 · 10 fracción kg/cm^3
Módulo de Young	210·105Ncm2
Coeficiente de Poisson	0.3
Límite elástico	24500Ncm2

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 fyint (N) y deflexión u_y (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 (t_n)

Flujo	i/tn	Error(%)	Max. Resid.	fyint	uy	Estado
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

Figura 12.

Campo de desplazamientos verticales en el arco elastoplástico. Fuerza en el nodo central igual a 6.6 · 10⁴ 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 · 10⁴ 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.

Figura 11.

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.

Figura 13.

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 × 10⁴. 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 14–22.

Figura 14.

Desplazamiento u_x (mm) obtenido con la técnica propuesta.

Figura 15.

Desplazamiento vertical u_y (mm) mediante Newton-Raphson y división de carga externa.

Figura 16.

Desplazamiento vertical u_y (mm) mediante la técnica propuesta.

Figura 17.

Tensión de cortadura σ_xy (N/mm²) mediante Newton-Raphson y división de carga externa.

Figura 18.

Tensión de cortadura σ_xy (N/mm²) mediante la técnica propuesta.

Figura 19.

Tensión σ_xx (N/mm²) mediante Newton-Raphson y división de carga externa.

Figura 20.

Tensión σ_xx (N/mm²) mediante la técnica propuesta.

Figura 21.

Tensión σ_yy (N/mm²) mediante Newton-Raphson y división de carga externa.

Figura 22.

Tensión σ_yy (N/mm²) 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í.

11

Conclusiones

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.

Agradecimientos

El primer autor está agradecido por el soporte recibido por Rockfield Software Ltd.