A stabilized mixed explicit formulation for plasticity with strain localization

Resumen

Este artículo presenta la aplicación de la formulación mixta estabilizada explícita en desplazamientos y deformaciones (MEX-FEM) [23,24] para la solución de problemas no lineales de la mecánica de sólidos con localización de deformaciones. A fin de emplear el mismo orden lineal de interpolación para el campo de los desplazamientos y deformaciones, nuestra formulación emplea el método de las sub-escalas variacionales. Comparada con la formulación estándar en desplazamientos, la formulación propuesta proporciona mejores campos de deformaciones y tensiones, y es capaz de abordar situaciones quasi-incompresibles. En este trabajo se investigan los efectos que tienen las deformaciones y tensiones mejoradas en los modelos de plasticidad de Mohr-Coulomb y Drucker Prager, incluyendo el fenómeno de la localización de las deformaciones. Los ejemplos numéricos validan la capacidad de la formulación propuesta para predecir correctamente los mecanismos de fallo, cargas últimas y la dirección de la banda de localización, virtualmente independientes de la malla utilizada y sin necesidad de emplear un algoritmo de rastreo.

Abstract

This paper presents the application of stabilized mixed explicit strain/displacement formulation (MEX-FEM) [23,24] for solving non-linear plasticity problems in solid mechanics with strain localization. In order to use the same linear interpolation order for displacements and strains, the formulation uses the variational subscales method. Compared to the standard irreducible formulation, the proposed formulation yields improved strain and stress fields, and it is capable of addressing nearly incompressible situations. This work investigates the effects of the improved strain and stress fields in problems involving strain softening and localization leading to failure for the Mohr-Coulomb and Drucker Prager plasticity models. Numerical examples validate the ability of the proposed formulation to correctly predict failure mechanisms with localized patterns of strain, virtually free of mesh dependence and without using tracking algorithm.

Palabras clave

Formulación mixta explícita;Estabilización;Plasticidad compresible e incompresible;Localización de deformaciones

Keywords

Explicit mixed formulation;Stabilization;Compressible and incompressible plasticity;Strain localization

1. Introducción

En trabajos anteriores [23] ;  [24] los autores han propuesto una formulación mixta explícita de elementos finitos (MEX-FEM) en desplazamientos y deformaciones para abordar problemas cuasi-estáticos y dinámicos en elasticidad y plasticidad compresible y cuasi-incompresible. La formulación propuesta usa elementos con independiente e igual orden lineal de interpolación para las variables involucradas estabilizada mediante sub-escalas variacionales (VMS)[20]. El objetivo de esta formulación es mejorar la precisión del campo de deformaciones y tensiones del problema discreto. Esta precisión es crucial para abordar problemas de plasticidad que presentan localización de deformaciones y ablandamiento del material [6]. En muchos casos de interés, con flujo plástico fundamentalmente isocórico, el comportamiento del modelo elasto-plástico tiende a la incompresibilidad en la medida en que las deformaciones plásticas predominan sobre las deformaciones elásticas, lo que lleva a un bloqueo volumétrico local en la vecindad de banda de localización. Los elementos finitos basados en la formulación irreducible, sobre todo los elementos de bajo orden, presentan dificultades para representar este comportamiento. Esto se traduce en la mayoría de los casos en una predicción errónea de la carga última y resultados dependiente de la malla.

Las formulaciones mixtas en desplazamientos-presión (u − p) son una alternativa viable y atractiva para solventar estos problemas. Cervera et al [4]; [5]; [7]; [17] ;  [33] han empleado esta formulación para la solución de problemas de plasticidad J2 con ablandamiento obteniendo resultados satisfactorios y virtualmente libres de la dependencia de la malla. Sin embargo, como las deformaciones desviadoras son calculadas por diferenciación del campo de los desplazamientos, se obtiene la misma velocidad de convergencia para las deformaciones que en el caso irreducible. Para mejorar el campo de las deformaciones desviadoras, Cervera et al [8] ;  [9] han introducido una formulación mixta estabilizada implícita en desplazamientos y deformaciones (u − ɛ) y seguidamente la han aplicado para abordar problemas de localización en el modelo de plasticidad de Drucker-Prager [6].

El objetivo de este trabajo es extender la aplicabilidad de MEX-FEM al rango elasto-plástico y demostrar que una formulación mixta estabilizada explícita es también una opción viable para resolver satisfactoriamente problemas cuasi-estáticos no lineales que involucran localización de las deformaciones. La organización del presente artículo es la siguiente. En la sección 2 se presentan el planteamiento mixto estándar en plasticidad y la extensión de la formulación propuesta al rango elasto-plástico. En la sección 3 se describen brevemente los modelos clásicos de plasticidad de Mohr-Coulomb y Drucker-Prager. Se detallan la integración constitutiva, el algoritmo de retorno y algunos aspectos generales sobre la dirección de la banda de localización. Por último, la sección 4 presenta una serie de ejemplos numéricos que validan el buen comportamiento de los elementos de la formulación propuesta.

2. Formulación mixta estabilizada explícita en desplazamiento-deformación u − ɛ en elasto-plasticidad

2.1. Formulación mixta u − ɛ en plasticidad

En problemas no lineales de la mecánica de sólidos, las deformaciones ɛ pueden considerarse como variables independientes adicionales al campo de los desplazamientos u (y sus derivadas temporales). En este caso, la forma fuerte del problema continuo se escribe como: dados los valores prescritos de las tracciones externas ${\textstyle {\overline {\boldsymbol {t}}}={\boldsymbol {\sigma }}{\boldsymbol {n}}:\partial {\Omega }_{\sigma }\rightarrow {\mathbb {R} }^{n_{dim}}}$, los desplazamientos u y aceleraciones ${\textstyle {\ddot {\boldsymbol {u}}}}$ en ∂Ω_u y las fuerzas por unidad de volumen ${\textstyle {\boldsymbol {b}}:\Omega \rightarrow {\mathbb {R} }^{n_{dim}}}$, hallar el campo de los desplazamientos, velocidades, aceleraciones y deformaciones en cualquier instante de tiempo ${\textstyle t\in I}$ siendo ${\textstyle I=(0,T)}$ el intervalo de tiempo de interés, tales que:

$\nabla \cdot {\boldsymbol {\sigma }}+{\boldsymbol {b}}=\rho {\ddot {\boldsymbol {u}}}\quad en\quad \Omega$

(1)

donde ${\textstyle \Omega \subset {\mathbb {R} }^{n_{dim}}}$ es el volumen ocupado por el sólido en un espacio de n_dim dimensiones; ρ denota la densidad del material y σ(ɛ) es el tensor de tensiones de Cauchy.

En la teoría de plasticidad infinitesimal, las deformaciones totales ɛ pueden descomponerse en la suma de dos contribuciones, una parte elástica ɛ_e y otra plástica ɛ_p,

${\boldsymbol {\varepsilon }}={\boldsymbol {\varepsilon }}_{e}+$${\boldsymbol {\varepsilon }}_{p}$

(3)

La ecuación constitutiva en plasticidad es:

${\boldsymbol {\sigma }}={\boldsymbol {C}}_{0}:\left({\boldsymbol {\varepsilon }}-\right.$$\left.{\boldsymbol {\varepsilon }}_{p}\right)={\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}_{e}$

(4)

donde C₀ es el tensor constitutivo elástico inicial. El problema queda definido formulando una ley de evolución del tipo ${\textstyle {\dot {{\boldsymbol {\varepsilon }}_{p}}}={\dot {\boldsymbol {\varepsilon }}}_{p}({\boldsymbol {\sigma }})}$ para las deformaciones plásticas.

La ecuación (4) puede expresarse en función de las deformaciones totales como [6]:

${\boldsymbol {\sigma }}=\left[{\boldsymbol {C}}_{0}-{\frac {{\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}_{p}\otimes {\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}_{p}}{{\boldsymbol {\varepsilon }}:{\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}_{p}}}\right]:{\boldsymbol {\varepsilon }}=$${\boldsymbol {C}}({\boldsymbol {\varepsilon }},{\boldsymbol {\varepsilon }}_{p}):{\boldsymbol {\varepsilon }}$

(5)

siendo C(ɛ, ɛ_p) el tensor constitutivo secante no lineal. La forma secante, no habitual en la teoría de la plasticidad, se adopta para usar el formalismo de las referencias [6] y [23]. Esto no afecta a los algoritmos de integración de las deformaciones plásticas y las tensiones. Esta forma de definir el tensor secante elasto-plástico tiene las simetrías deseadas de un tensor constitutivo secante.

Además de las ecuaciones (1) y (2), las variables del problema u y ɛ están sujetas a unas condiciones iniciales en t = t₀, es decir; ${\textstyle {\boldsymbol {u}}\vert _{t=t_{0}}={\boldsymbol {u}}_{0}}$, ${\textstyle {\dot {\boldsymbol {u}}}\vert _{t=t_{0}}={\boldsymbol {v}}_{0}}$ y ${\textstyle {\boldsymbol {\varepsilon }}\vert _{t=t_{0}}={\boldsymbol {\varepsilon }}_{0}=}$${\nabla }^{s}{\boldsymbol {u}}_{0}$. Después de multiplicar las ecuaciones (1) y (2) por sus respectivas funciones de peso e integrando por partes la ecuación (1), la formulación variacional del problema es:

${\int }_{\Omega }{\nabla }^{s}{\boldsymbol {\omega }}:{\boldsymbol {\sigma }}\quad d\Omega +$${\int }_{\Omega }{\boldsymbol {\omega }}\cdot \rho {\ddot {\boldsymbol {u}}}\quad d\Omega =$${\int }_{\partial \Omega }{\boldsymbol {\omega }}\cdot {\overline {\boldsymbol {t}}}\quad d\Gamma +$${\int }_{\Omega }{\boldsymbol {\omega }}\cdot {\boldsymbol {b}}\quad d\Omega \quad \quad \forall {\boldsymbol {\omega }}$

(6)

donde ${\textstyle {\boldsymbol {\omega }}\in V}$ y ${\textstyle {\boldsymbol {\gamma }}\in T}$ son las funciones de prueba para el campo de los desplazamientos y el campo de las deformaciones, respectivamente; ${\textstyle V}$ y ${\textstyle T}$ son los espacios de los desplazamientos y deformaciones admisibles.

Para la definición discreta del problema, se introducen en las ecuaciones (6) y (7) un campo de desplazamientos discretos, u = u_h, aceleraciones discretas ${\textstyle {\ddot {\boldsymbol {u}}}={\ddot {\boldsymbol {u}}}_{h}}$ y deformaciones discretas ɛ = ɛ_h tales que:

${\int }_{\Omega }{\nabla }^{s}{\boldsymbol {\omega }}_{h}:{\boldsymbol {\sigma }}\quad d\Omega +$${\int }_{\Omega }{\boldsymbol {\omega }}_{h}\cdot \rho {\ddot {\boldsymbol {u}}}_{h}\quad d\Omega =$${\int }_{\partial \Omega }{\boldsymbol {\omega }}_{h}\cdot {\overline {\boldsymbol {t}}}\quad d\Gamma +$${\int }_{\Omega }{\boldsymbol {\omega }}_{h}\cdot {\boldsymbol {b}}\quad d\Omega \quad \quad \forall {\boldsymbol {\omega }}_{h}$

(8)

en donde ${\textstyle {\boldsymbol {\omega }}_{h}\in V_{h}}$ y ${\textstyle {\boldsymbol {\gamma }}_{h}\in T_{h}}$ son las funciones de interpolación definidas sobre el espacio de los elementos finitos ${\textstyle V_{h}}$ y ${\textstyle T_{h}}$, respectivamente. Es de interés en este trabajo el empleo de funciones de interpolación lineales tanto para el campo de los desplazamientos como para el campo de las deformaciones. Sin embargo, el empleo de las mismas funciones de interpolación en ambos campos no cumple la condición inf-sup LBB [2] ;  [18], por lo que la solución obtenida con las ecuaciones (8) y (9) es inestable y presenta oscilaciones en el campo de los desplazamientos. A fin evitar la condición LBB, en este trabajo se emplea el método de las Multi-Escalas Variacionales (VMS)[21].

2.2. Formulación mixta estabilizada explícita (MEX-FEM) en elasto-plasticidad

En el método VMS se introducen dos niveles de resolución (escalas), una que se logra captar mediante los elementos finitos (escala gruesa), y otra fina, que corresponde a la parte que la malla no logra captar (sub-escala). La solución continua del problema contiene entonces componentes de ambas escalas. Por lo tanto, empleando VMS, los desplazamientos y las deformaciones se aproximan como:

${\boldsymbol {u}}={\boldsymbol {u}}_{h}+{\tilde {\boldsymbol {u}}}$

(10)

Asimismo, derivando en el tiempo,

${\dot {\boldsymbol {u}}}={\dot {\boldsymbol {u}}}_{h}+{\dot {\tilde {\boldsymbol {u}}}}$

(12)

donde ${\textstyle {\boldsymbol {\varepsilon }}_{h}\in T_{h}}$ y ${\textstyle ({\boldsymbol {u}}_{h},{\dot {\boldsymbol {u}}}_{h},{\ddot {\boldsymbol {u}}}_{h})\in V_{h}}$ son las deformaciones, desplazamientos, velocidades y aceleraciones en la escala gruesa; ${\textstyle {\tilde {\boldsymbol {\varepsilon }}}\in {\tilde {T}}}$ y ${\textstyle ({\tilde {\boldsymbol {u}}},{\dot {\tilde {\boldsymbol {u}}}},{\ddot {\tilde {\boldsymbol {u}}}})\in {\tilde {V}}}$ son las deformaciones y desplazamientos (con sus derivadas temporales) de la sub-escala y ${\textstyle (T_{h},V_{h})}$ e ${\textstyle ({\tilde {T}},{\tilde {V}})}$ son los espacios de las funciones admisibles de las deformaciones y los desplazamientos para la escala gruesa y las sub-escalas, respectivamente. Para el caso elástico lineal, las tensiones pueden escribirse como:

${\boldsymbol {\sigma }}={\boldsymbol {C}}_{0}:\left({\boldsymbol {\varepsilon }}_{h}+\right.$$\left.{\tilde {\boldsymbol {\varepsilon }}}\right)={\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}_{h}+$${\boldsymbol {C}}_{0}:{\tilde {\boldsymbol {\varepsilon }}}=$${\boldsymbol {\sigma }}_{h}+{\tilde {\boldsymbol {\sigma }}}$

(14)

Introduciendo las ecuaciones (10), (12), (13), (11) y (14) en las ecuaciones (6) y (7) se obtiene:

${\int }_{\Omega }{\nabla }^{s}{\boldsymbol {\omega }}_{h}:{\boldsymbol {\sigma }}_{h}\quad d\Omega +$${\int }_{\Omega }{\nabla }^{s}{\boldsymbol {\omega }}_{h}:{\tilde {\boldsymbol {\sigma }}}\quad d\Omega +$${\int }_{\Omega }{\boldsymbol {\omega }}_{h}\cdot \rho \left({\ddot {\boldsymbol {u}}}_{h}+\right.$$\left.{\ddot {\tilde {\boldsymbol {u}}}}\right)\quad d\Omega =$$F_{ext}({\boldsymbol {\omega }}_{h})\quad \forall {\boldsymbol {\omega }}_{h}$

(15)

donde ${\textstyle F_{ext}({\boldsymbol {\omega }}_{h})={\int }_{\partial \Omega }{\boldsymbol {\omega }}_{h}\cdot {\overline {\boldsymbol {t}}}\quad d\Gamma +}$${\int }_{\Omega }{\boldsymbol {\omega }}_{h}\cdot {\boldsymbol {b}}\quad d\Omega$. Si bien las ecuaciones (15) y (17) se plantean y se resuelven en el espacio de los elementos finitos, no así las ecuaciones (16) y (18), por lo que se tiene que recurrir a un procedimiento aproximado a fin evaluar los términos de las sub-escalas en las ecuaciones (15) y (17). Aquí se emplea el método de las Sub-escalas Ortogonales OSS.

Las ecuaciones (15) y (16) requieren un método de integración temporal. En este trabajo se usa el método de integración explícita de Diferencias Centradas (CD) en ambas escalas. Los detalles de MEX-FEM para el caso elástico compresible y cuasi-incompresible se dan en la referencia [24].

La extensión de MEX-FEM al rango elasto-plástico radica fundamentalmente en dos aspectos. En primer lugar, el tensor de tensiones σ es no lineal en ɛ. En segundo lugar, en el rango elasto-plástico los parámetros de estabilización varían en la medida en que se desarrolla el flujo plástico.

Suponiendo que la contribución de la sub-escala es pequeña respecto a la escala resoluble, es decir ${\textstyle \parallel {\tilde {\boldsymbol {\varepsilon }}}\parallel <<\parallel {\boldsymbol {\varepsilon }}_{h}\parallel }$ (véase [6]; [8]; [23] ;  [24]), es posible aproximar

${\boldsymbol {\varepsilon }}_{p}={\boldsymbol {\varepsilon }}_{p}({\boldsymbol {\sigma }})\simeq {\boldsymbol {\varepsilon }}_{p}({\boldsymbol {\sigma }}_{h})$

(19)

lo que implica que las tensiones continuas son:

${\boldsymbol {\sigma }}={\boldsymbol {C}}_{0}:\left({\boldsymbol {\varepsilon }}_{h}+\right.$$\left.{\tilde {\boldsymbol {\varepsilon }}}-{\boldsymbol {\varepsilon }}_{p}({\boldsymbol {\sigma }}_{h})\right)=$${\boldsymbol {\sigma }}_{h}+{\tilde {\boldsymbol {\sigma }}}$

(20)

En el método VMS, las sub-escalas no están definidas de forma local, sino de forma variacional. Por tanto, para utilizar la sub-escala de deformaciones en el cálculo de la ecuación constitutiva, y en particular en el cálculo de las deformaciones plásticas, hay dos opciones: a) localizar su valor por algún procedimiento heurístico o b) despreciar su contribución. Ambas variantes producen resultados comparables lo que valida el hecho de la contribución de la sub-escala de las deformaciones es “pequeña”. Obsérvese que la expresión (20) es análoga a la expresión (14). De esta manera se tiene el mismo formalismo y procedimiento que en el caso elástico [24].

Para calcular el valor de las sub-escalas se escoge el espacio ortogonal al espacio de los elementos finitos (Método de las Sub-escalas Ortogonales (OSS) [12]; [13] ;  [14]. Por practicidad se introduce el operador de proyección ortogonal ${\textstyle P_{h}^{\perp }(\bullet )=(\bullet )-P_{h}(\bullet )}$. Dado que la sub-escala de los desplazamientos es dinámica, se integra ésta en el tiempo usando el método explícito de CD   e introduciendo un coeficiente de disipación ${\textstyle {\tilde {\xi }}\in [0,1]}$ se tiene que:

${\tilde {\boldsymbol {u}}}^{n+1}\simeq {\tau }_{u_{t}}{\frac {\rho }{\Delta t^{2}}}\left((2-\right.$$\left.{\tilde {\xi }}){\tilde {\boldsymbol {u}}}^{n}-(1-\right.$$\left.{\tilde {\xi }}){\tilde {\boldsymbol {u}}}^{n-1}\right)+$${\tau }_{u_{t}}P^{\perp }\left(\nabla \cdot {\boldsymbol {\sigma }}_{h}^{n}\right)$

(21)

Los términos ${\textstyle {\tau }_{u_{t}}}$ y τ_ɛ son los parámetros de estabilización dados por:

${\tau }_{u_{t}}={\left({\frac {\rho }{\Delta t^{2}}}+B(\mu ){\frac {{\mu }_{0}}{c_{u}h_{e}^{2}}}\right)}^{-1}\quad y\quad {\tau }_{\varepsilon }=$$c_{\varepsilon }B(\mu ){\frac {h_{e}}{L_{0}}}$

(23)

donde c_u > 0 y c_ɛ > 0 son constantes algorítmicas adimensionales, h_e es la longitud del elemento finito, L₀ es la longitud característica del problema, μ₀ = 2G es el módulo inicial de rigidez al corte y β(μ) = μ/μ₀, siendo ${\textstyle \mu ={\frac {\parallel {\boldsymbol {S}}_{h}\parallel }{\parallel dev({\boldsymbol {\varepsilon }}_{h})\parallel }}}$ el módulo de corte secante [28] ;  [33].

El coeficiente de corte secante decrece rápidamente a medida que localizan las deformaciones, lo cual puede causar oscilaciones es espurias en la sub-escala de los desplazamientos. Una manera de sortear este problema es reemplazar el coeficiente de corte efectivo μ   por ${\textstyle {\overline {\mu }}=(1-{\overline {\eta }}){\mu }_{old}+{\overline {\eta }}{\mu }_{new}}$ donde ${\textstyle {\overline {\eta }}\in \lbrace 0,1\rbrace }$ puede interpretarse como un coeficiente de retardo.

Cuando se desarrolla el flujo plástico y se produce la localización de deformaciones, las deformaciones desviadoras totales y, particularmente, las deformaciones desviadoras plásticas, predominan sobre las deformaciones volumétricas y se tiene un comportamiento cuasi-incompresible. La razón fundamental para introducir la subescala de los desplazamientos es poder estabilizar la tensión media ${\textstyle p={\frac {1}{3}}tr({\boldsymbol {\sigma }})}$ en situaciones cuasi-incompresibles. El término que estabiliza el campo de la tensión media es el gradiente ∇p_h contenido dentro de ∇ · σ_h = ∇ · S_h + ∇ p_h siendo S_h el tensor desviador de tensiones. Se ha demostrado efectivo sustituir el último término de la ecuación (21), ${\textstyle P_{h}^{\perp }(\nabla \cdot {\boldsymbol {\sigma }}_{h})}$, por ${\textstyle P_{h}^{\perp }(\nabla p_{h})}$, quedando

${\tilde {\boldsymbol {u}}}^{n+1}\simeq {\tau }_{u_{t}}{\frac {\rho }{\Delta t^{2}}}\left((2-\right.$$\left.{\tilde {\xi }}){\tilde {\boldsymbol {u}}}^{n}-(1-\right.$$\left.{\tilde {\xi }}){\tilde {\boldsymbol {u}}}^{n-1}\right)+$${\tau }_{u_{t}}\left(\nabla p_{h}^{n}-P_{h}\left(\nabla p_{h}^{n}\right)\right)$

(24)

2.3. Formulación matricial de la formulación mixta explícita

Las funciones de prueba ω_h para el espacio de los desplazamientos y γ_h para el espacio de las deformaciones, en las ecuaciones (15) y (17), respectivamente, constituyen las funciones de interpolación. Para el caso de un tetraedro lineal de cuatro nodos, empleando las las mismas funciones lineales de interpolación para los desplazamientos y deformaciones se tiene que [24]:

${\boldsymbol {u}}_{h}=\sum _{i=1}^{4}{\boldsymbol {N}}_{u}^{i}{\boldsymbol {u}}_{h}^{i}\quad y\quad {\boldsymbol {\varepsilon }}_{h}=$$\sum _{i=1}^{4}{\boldsymbol {N}}_{\boldsymbol {\varepsilon }}^{i}{\boldsymbol {\varepsilon }}_{h}^{i}$

(25)

siendo ${\textstyle {\boldsymbol {N}}_{u}^{i}}$ la submatriz diagonal de 3 × 3 ${\textstyle ({\boldsymbol {N}}_{u}^{i}=diag\lbrace N_{i},N_{i},N_{i}\rbrace )}$ de las funciones de interpolación del campo de los desplazamientos en el nodo i  ${\textstyle ({\boldsymbol {u}}_{h}^{i}=\lbrace u_{x,h}^{i},u_{y,h}^{i},u_{z,h}^{i}\rbrace )}$ y ${\textstyle {\boldsymbol {N}}_{\boldsymbol {\varepsilon }}^{i}}$ la submatriz diagonal de 6 × 6 ${\textstyle ({\boldsymbol {N}}_{\boldsymbol {\varepsilon }}^{i}=diag\lbrace N_{i},N_{i},N_{i},N_{i},N_{i},N_{i}\rbrace )}$ de las funciones de interpolación del campo de las deformaciones en el nodo i  ${\textstyle ({\boldsymbol {\varepsilon }}_{h}^{i}=\lbrace {\boldsymbol {\varepsilon }}_{xx,h}^{i},{\boldsymbol {\varepsilon }}_{yy,h}^{i},{\boldsymbol {\varepsilon }}_{zz,h}^{i},{\boldsymbol {\varepsilon }}_{xy,h}^{i},{\boldsymbol {\varepsilon }}_{yz,h}^{i},{\boldsymbol {\varepsilon }}_{xz,h}^{i}\rbrace )}$. Por otro lado, en la ecuación (15) aparece el término ${\textstyle {\nabla }^{s}{\boldsymbol {\omega }}_{h}={\boldsymbol {B}}_{u}^{T}}$ que no es más que el clásico operador matricial de gradiente simétrico discreto, en el cual B_u = [B₁B₂B₃B₄], donde cada sub-matriz B_i se escribe como:

${\boldsymbol {B}}_{i}={\left({\begin{array}{cccccc}{\frac {\partial N_{i}}{\partial X}}&0&0&{\frac {\partial N_{i}}{\partial Y}}&0&{\frac {\partial N_{i}}{\partial Z}}\\0&{\frac {\partial N_{i}}{\partial Y}}&0&{\frac {\partial N_{i}}{\partial X}}&{\frac {\partial N_{i}}{\partial Z}}&0\\0&0&{\frac {\partial N_{i}}{\partial Z}}&0&{\frac {\partial N_{i}}{\partial Y}}&{\frac {\partial N_{i}}{\partial X}}\\&&&&&\end{array}}\right)}^{T}$

(26)

Asimismo, el operador de divergencia ∇ · γ_h en la ecuación (17) toma la forma de ${\textstyle {\boldsymbol {B}}_{\varepsilon }=[{\boldsymbol {B}}_{1}^{T}{\boldsymbol {B}}_{2}^{T}{\boldsymbol {B}}_{3}^{T}{\boldsymbol {B}}_{4}^{T}]}$.

El término inercial de la ecuación (15) da lugar al vector de fuerzas de inercia, ${\textstyle {\boldsymbol {M}}{\ddot {\boldsymbol {u}}}_{h}}$, donde ${\textstyle {\boldsymbol {M}}={\int }_{{\Omega }_{e}}{\boldsymbol {N}}_{u}^{T}\rho {\boldsymbol {N}}_{u}d\Omega }$ es la matriz de masa del sistema. Se introduce el amortiguamiento estructural usando el clásico amortiguamiento de Rayleigh proporcional a la matriz de masa, es decir, D = α''M. La integración temporal explícita es efectiva si las matrices de masa M y de amortiguamiento D son diagonales. En tal caso, la forma matricial mixta explícita se escribe como [23] ;  [24]:

${\boldsymbol {\varepsilon }}_{h}^{n}=P_{h}\left({\nabla }^{s}{\boldsymbol {u}}^{n}\right)-$$P_{h}\left(\left(\nabla \cdot {\boldsymbol {\gamma }}_{h}\right)\cdot {\tilde {\boldsymbol {u}}}^{n}\right)$

(27)

siendo ${\textstyle {\boldsymbol {F}}_{int}^{stab,n}={\int }_{\Omega }{\boldsymbol {B}}_{\boldsymbol {u}}^{T}{\boldsymbol {\sigma }}({\boldsymbol {\varepsilon }}_{n}^{stab})\quad d\Omega }$ y ${\textstyle {\boldsymbol {F}}_{ext}^{n}={\int }_{\Gamma }{\boldsymbol {N}}_{\boldsymbol {u}}{\overline {\boldsymbol {t}}}d\Gamma +}$${\int }_{\Omega }{\boldsymbol {N}}_{\boldsymbol {u}}{\boldsymbol {b}}d\Omega$ los vectores de fuerzas internas estabilizadas y fuerzas externas, respectivamente. Las aceleraciones y las velocidades ${\textstyle [{\ddot {\boldsymbol {u}}}_{h}^{n},{\dot {\boldsymbol {u}}}_{h}^{n}]}$ se aproximan mediante diferencias centradas en la forma:

${\ddot {\boldsymbol {u}}}_{h}^{n}={\frac {1}{\Delta t^{2}}}\left({\boldsymbol {u}}_{h}^{n+1}-\right.$$\left.2{\boldsymbol {u}}_{h}^{n}+{\boldsymbol {u}}_{h}^{n-1}\right)$

(30)

El algortimo de resolución de las ecuaciones mixta explícita sigue los pasos siguientes [24]:

• 1. Cómputo de las fuerzas internas estabilizadas ${\textstyle {\boldsymbol {F}}_{int}^{stab,n}}$.
• 2. Cómputo de los desplazamientos en t = tⁿ⁺¹:

${\boldsymbol {\mbox{u}}}^{n+1}={\left[2{\boldsymbol {M}}+\Delta t{\boldsymbol {D}}\right]}^{-1}[(4{\boldsymbol {M}}{\boldsymbol {\mbox{u}}}^{n}-$$[2{\boldsymbol {M}}-\Delta t{\boldsymbol {D}}]{\boldsymbol {\mbox{u}}}^{n-1}+$$2\Delta t^{2}({\boldsymbol {F}}_{ext}^{n}-{\boldsymbol {F}}_{int}^{stab,n})]$

• 3. Evaluación de las sub-escalas de los desplazamientos ${\textstyle {\tilde {\boldsymbol {u}}}^{n+1}}$ empleando la ecuación (24).
• 4. Actualización de las deformaciones ${\textstyle {\boldsymbol {\varepsilon }}_{h}^{n+1}}$ con la ecuación (27) en tⁿ⁺¹.
• 5. Ir a siguiente paso. La ventaja de emplear un esquema mixto explícito de integración es que no requiere la resolución iterativa de un sistema no diagonal de ecuaciones no lineales, ni el cómputo del tensor constitutivo secante no lineal C. Sin embargo, el esquema propuesto es condicionalmente estable, por lo que el paso de tiempo Δt está sujeto al límite de Courant. No obstante, se ha demostrado numéricamente que el paso de tiempo en la formulación mixta explícita es mayor que el de la formulación irreducible explícita [23].

3. Modelos clásicos de plasticidad

3.1. Modelo de plasticidad de Mohr-Coulomb y Drucker-Prager

El criterio de Mohr-Coulomb (MC) se emplea para describir el fallo en materiales friccionales y geo-materiales en general. El comportamiento de estos materiales se caracterizan por la dependencia de la cohesión efectiva con la presión. Las deformaciones plásticas son el resultado del deslizamiento relativo y la fricción entre las partículas. De acuerdo con este criterio, el flujo plástico comienza cuando cierta combinación del esfuerzo cortante τ y el esfuerzo normal σ_n alcanzan un valor crítico:

$\vert \tau \vert =c-{\sigma }_{n}tan\phi$

(32)

siendo c = c(ɛ_p) ≥ 0 la cohesión y 0 ≤ ϕ ≤ π/2 el ángulo de fricción interna. En el caso de que ϕ = 0 se tiene el modelo de plasticidad de Tresca. La localización de deformaciones se producen si la cohesión descrece a la medida en que se incrementan las deformaciones plásticas. Entonces:

$\Phi \left(\tau ,{\sigma }_{n},c({\boldsymbol {\varepsilon }}_{p}),\phi \right)=$$\vert \tau \vert -c({\boldsymbol {\varepsilon }}_{p})+$${\sigma }_{n}tan\phi =0$

(33)

El criterio de MC también puede expresarse en función de los invariantes del tensor de tensiones o más habitualmente en 6 funciones expresadas en el espacio de Haigh-Westergaard (HW) (o de tensiones principales) que combinadas forman la representación en superficie múltiples del modelo de MC. La superficie principal está dada por:

${\Phi }_{1}\left({\sigma }_{1},{\sigma }_{3},c({\boldsymbol {\varepsilon }}_{p}),\phi \right)=$$({\sigma }_{1}-{\sigma }_{3})+({\sigma }_{1}+{\sigma }_{3})sin\phi -$$2c({\boldsymbol {\varepsilon }}_{p})cos\phi =0$

(34)

siendo σ₁ ≥ σ₂ ≥ σ₃ las tensiones principales. En el espacio de HW, el modelo de MC es una pirámide hexagonal cuyo ápice se localiza en ${\textstyle {\sqrt {3}}ccot\phi }$ tal como se muestra en la figura 1.

Figura 1. Modelo de Mohr-Coulomb.

El criterio de Drucker-Prager (1952) fue propuesto por Drucker y Prager como una aproximación suavizada   del criterio de MC. Es una modificación del modelo de Von Mises en la que se incluye la dependencia con la presión. Este modelo sugiere que el fallo en el material comienza cuando el segundo invariante del tensor desviador de tensiones ${\textstyle J_{2}={\frac {1}{2}}{\boldsymbol {S}}:{\boldsymbol {S}}}$ y la tensión hidrostática alcanzan una combinación crítica:

$\Phi \left({\boldsymbol {\sigma }},c({\boldsymbol {\varepsilon }}_{p}),\phi \right)=$${\sqrt {J_{2}({\boldsymbol {\sigma }})}}+\eta p-\varsigma c({\boldsymbol {\varepsilon }}_{p})=$$0$

(35)

El criterio de DP, mostrado en la figura 2, es en un cono cuyo eje coincide con el eje hidrostático. Nótese que para ϕ = 0 se recupera el modelo de plasticidad incompresible de Von Mises. Los parámetros η y ς se escogen de acuerdo a la aproximación de MC que se desea realizar. A fin de predecir idénticas cargas últimas para los modelos de MC y DP en un estado de deformación plana, η y ς se toman tales que [26]:

Figura 2. Modelo de Drucker-Prager.

$\eta ={\frac {3tan\phi }{\sqrt {9+12{tan}^{2}\phi }}}\quad y\quad \varsigma =$${\frac {3}{\sqrt {9+12{tan}^{2}\phi }}}$

(36)

El modelo de Drucker Prager (DP) es isocórico para ángulo de fricción nulo y cuasi-incompresible para ángulos de fricción pequeños. De igual manera, el modelo de Mohr-Coulomb (MC) es también poco compresible para ángulos de fricción pequeños. Por esta razón se introduce la sub-escala de los desplazamientos para tener un esquema numérico robusto y estable en estas condiciones.

3.2. Integración constitutiva, algoritmo de retorno y de ablandamiento

El objetivo de la integración constitutiva es encontrar los valores actuales de las variables plásticas. Para integrar numéricamente un modelo de plasticidad se emplean habitualmente algoritmos implícitos de retorno (return mapping). El cálculo de la velocidad de deformación plástica y la deformación plástica efectiva en plasticidad asociada con superficie múltiples se realiza mediante la regla de Koiter [31]:

${\dot {\boldsymbol {\varepsilon }}}_{p}=\sum _{\alpha =1}^{n_{act}}{\dot {\lambda }}_{\alpha }{\frac {\partial {\Phi }_{\alpha }}{\partial {\boldsymbol {\sigma }}}}=$$\sum _{\alpha =1}^{n_{act}}{\dot {\lambda }}_{\alpha }{\boldsymbol {m}}_{\alpha }$

(37)

donde ${\textstyle {\dot {\lambda }}_{\alpha }}$ son los parámetros de consistencia plástica y n_act son las superficies ${\textstyle J}$activas. El problema consiste en determinar ɛ_p, λ_α y las superficies involucradas ${\textstyle J:={\Phi }_{\alpha }\forall \alpha \in \lbrace 1,2..n_{act}\rbrace }$. Asimismo, los parámetros de consistencia plástica siguen las condiciones complementarias de Kuhn-Tucker extendidas a modelos de plasticidad con superficies múltiples:

${\dot {\lambda }}_{\alpha }\geq 0,\quad {\dot {\lambda }}_{\alpha }{\Phi }_{\alpha }\left({\boldsymbol {\sigma }},c({\boldsymbol {\varepsilon }}_{p})\right)\leq 0\quad y\quad {\dot {\lambda }}_{\alpha }{\dot {\Phi }}_{\alpha }\left({\boldsymbol {\sigma }},c({\boldsymbol {\varepsilon }}_{p})\right)=$$0$

(39)

Los procedimientos de la integración constitutiva del modelo de MC y DP empleados en este trabajo están descritos en la referencia [26].

Por otro lado, la energía disipada durante la formación de la banda de corte está intrínsecamente relacionada con la energía de fractura ${\textstyle G_{f}}$ por unidad de área. En un proceso uniaxial de carga, con ablandamiento y relajación total final de las tensiones, el trabajo plástico ${\textstyle W_{p}}$ realizado por unidad de volumen representa la energía dispada durante la carga plástica. Estas dos cantidades están relacionadas a través de la longitud característica del elemento finito h_e como:

$W_{p}={\frac {G_{f}}{h_{e}}}$

(40)

Suponiendo que la evolución de la cohesión sigue lo descrito en la tabla 1, el trabajo total en todo el proceso plástico incluyendo el completo ablandamiento, esto es, desde el estado elástico (t = 0, c = c₀ y ${\textstyle {\overline {\boldsymbol {\varepsilon }}}_{p}=0}$) hasta el desarrollo de deformaciones plásticas localizadas (t =∞, c = 0 y ${\textstyle {\overline {\boldsymbol {\varepsilon }}}_{p}\not =0}$) es igual a:

Tabla 1. Curva de ablandamiento de la cohesión

${\textstyle c={\begin{array}{ll}c_{0}\left(1-{\frac {H_{s}{\overline {\boldsymbol {\varepsilon }}}_{p}}{c_{0}}}\right)&\quad si\quad 0<{\overline {\boldsymbol {\varepsilon }}}_{p}<{\frac {c_{0}}{H_{s}}}\quad {\mbox{ablandamiento lineal}}\\0&\quad si\quad {\overline {\boldsymbol {\varepsilon }}}_{p}>{\frac {c_{0}}{H_{s}}}.\end{array}}}$
${\textstyle c=c_{0}exp\left({\frac {-2H_{s}}{c_{0}}}{\overline {\boldsymbol {\varepsilon }}}_{p}\right)\quad para\quad 0<{\overline {\boldsymbol {\varepsilon }}}_{p}<\infty \quad {\mbox{ablandamiento exponencial}}}$

$W_{p}={\int }_{t=0}^{t=\infty }{\dot {W}}_{p}dt={\int }_{t=0}^{t=\infty }\sigma :{\dot {\boldsymbol {\varepsilon }}}_{p}dt=$${\frac {c_{0}^{2}}{2H_{s}}}$

(41)

El escalar ${\textstyle H_{s}}$ mide la fragilidad del material y depende únicamente de las propiedades de los materiales y de la discretización empleada.

3.3. Orientación de la banda de localización. Generalidades

La localización de la deformación se manifiesta en los materiales elasto-plásticos como una banda de corte, una zona estrecha de intensas deformaciones a través del cual los campos de deformaciones son discontinuos. Varios autores [1]; [3]; [11]; [22]; [25]; [27]; [29]; [30] ;  [32] han encontrado soluciones analíticas y geométricas para la orientación de las bandas de discontinuidad ${\textstyle S}$ (véase fig. 3) en modelos elasto-plásticos empleando diferentes estrategias. Todos ellos basan sus soluciones en la llamada condición de localización, que implica la pérdida de elipticidad de la ecuación de balance en el caso estático o de hiperbolicidad en el caso dinámico, demostrando que es una condición necesaria para la aparición de discontinuidades débiles y posteriormente el fallo localizado. Por otra parte, Cervera et al [7]; [10] ;  [34] proponen una metodología distinta para encontrar el valor analítico de la orientación de banda de localización y es la empleada en este trabajo. Este procedimiento formula las condiciones de acotabilidad de las tensiones y decohesión total, que combinados, son las condiciones necesarias para la formación de la banda de corte. Según estas condiciones, la banda de localización en un modelo de plasticidad asociada no depende de las constantes elásticas, sino únicamente del estado tensional y del vector de flujo plástico. La tabla 2 muestra los valores analíticos del ángulo de localización θ_loc para los modelos de DP y MC en tensión y deformación plana. Cervera et al [6] han verificado numéricamente por EF los resultados analíticos de localización en un modelo de plasticidad J2 y DP empleando una formulación mixta implícita.

Figura 3. Ángulo de localización.

Tabla 2. Ángulo de localización teórico para los modelos de DP y MC en tensión y deformación plana

Modelo	Tensión Plana	Deformación Plana
DP	${\textstyle {tan}^{2}{\theta }_{loc}={\frac {{\sigma }_{1}-2{\sigma }_{2}-2J_{2}^{1/2}\eta }{2{\sigma }_{1}-{\sigma }_{2}+2J_{2}^{1/2}\eta }}}$	${\textstyle {tan}^{2}{\theta }_{loc}={\frac {{\sigma }_{1}-{\sigma }_{2}-2J_{2}^{1/2}\eta }{{\sigma }_{1}-{\sigma }_{2}+2J_{2}^{1/2}\eta }}}$
MC	${\textstyle {\theta }_{loc}={\begin{array}{ll}0^{0}&si{\mbox{flujo normal al plano}}\quad A-B.\\\pm \left({45}^{o}-{\frac {\phi }{2}}\right)&{\mbox{si flujo normal al plano}}\quad A-C.\end{array}}}$	${\textstyle {\theta }_{loc}=\pm \left({45}^{o}-{\frac {\phi }{2}}\right)}$

La formulación mixta considera la ecuación geométrica (7) en forma débil, mientras que la formulación irreducible lo hace en forma fuerte. Por eso las deformaciones nodales son grados de libertad independientes de los desplazamientos. Esta mejora de la aproximación de la cinemática discreta respecto de la formulación irreducible permite la convergencia local de las deformaciones en situaciones cuasi-singulares, evitando la dependencia de orientación de la malla. Esta es la razón por la que la formulación mixta funciona en problemas de localización. Por otro lado, su velocidad de convergencia en deformaciones es siempre superior a la de las formulaciones irreducibles. Por tanto, se puede esperar que sus resultados sean mejores para cualquier grado de refinamiento.

Hay que señalar que las sub-escalas ortogonales son necesarias por las razones de estabilidad, al elegir interpolaciones lineales para despazamientos y deformaciones que no cumplen la condición inf-sup. En caso de que se usen interpolaciones mixtas compatibles no es necesaria la estabilización.

4. Simulaciones numéricas

La formulación MEX-FEM se aplica a continuación en una serie de simulaciones numéricas de problemas de plasticidad compresible e incompresible. Los análisis se realizan empleando elementos triangulares y tetraédricos con igual orden lineal interpolación para los desplazamientos y deformaciones. El tamaño del elemento finito se calcula como ${\textstyle h_{e}={\left(4/\pi \cdot A_{e}\right)}^{\frac {1}{2}}}$ siendo A_e el área del elemento finito correspondiente. En 3D   se toma ${\textstyle h_{e}={\left(6/\pi \cdot V_{e}\right)}^{\frac {1}{3}}}$ donde V_e el volumen del elemento finito. Se asume un ablandamiento exponencial para la cohesión. A fin de obtener una respuesta cuasi-estática en el tiempo se emplea un amortiguamiento proporcional con relajación dinámica. Asimismo, para evaluar los parámetros de estabilización se toma c_u = [0.25, 5] y c_ɛ = 1.0. Los análisis se han realizado con el programa de elementos finitos KRATOS [15] ;  [16], desarrollado en el Centro Internacional de Métodos Numéricos (CIMNE). Como pre y post-procesador se ha utilizado GiD [19], también desarrollado en CIMNE.

4.1. Zapata de Prantl. Plasticidad perfecta

Este ejemplo bidimensional en deformación plana se usa a menudo para verificar mecanismos de colapso y cargas últimas en modelos de plasticidad. Las figuras 4 muestran la geometría y la malla empleada en el análisis. La propiedad del material son: densidad ρ = 10⁴ Kg/m³, módulo de Young E = 10⁷ KPa, coeficiente de Poisson ν = 0.48, cohesión inicial c₀ = 490 KPa y ángulo de fricción ϕ = 20^o. Este ejemplo se analiza empleando el modelo de plasticidad perfecta de MC y DP con η y ς dados en la ecuación (36). Dada la simetría del problema, sólo se analiza la mitad del dominio (la mitad derecha). Las dimensiones L y B y longitud característica del problema L₀ se toman como 5m y 1m, respectivamente 5m. El dominio se discretiza empleando una malla no estructurada de 2438 nodos y 4731 elementos lineales, tanto para la formulación irreducible (FI − P1), como para la presente formulación (MEX − FEMP1P1). Se impone una velocidad instantánea y constante en el tiempo de 10⁻³m/s. Asimismo, el problema se analiza estáticamente empleando los elementos de orden cuadrático de la formulación irreducible (FI − P2). La solución analítica es P_lim/c₀ = 14.8, donde P_lim es la fuerza total (reacción) aplicada.

Figura 4. Geometría y malla empleada en el ensayo de Prantl.

Las figuras 5a y 5b muestran la distribución de los campos de presión y de deformación plástica equivalente obtenidos con la formulación propuesta en el modelo de MC. Se obtuvieron distribuciones similares con el modelo de DP. Nótese que el campo de las presiones está estabilizado. El mecanismo de colapso predicho por esta formulación es conciso sin ninguna ramificación espuria. La banda de localización se forma de acuerdo a la solución clásica virtualmente libre de la dependencia de la malla. Las figuras 5c y 5d muestran las distribuciones obtenidas con la formulación irreducible (FI − P1). Obsérvese que ésta no produce soluciones satisfactorias, presentándose una ramificación espuria en la banda de localización y oscilaciones locales en el campo de las presiones. El elemento cuadrático de la formulación irreducible (FI − P2) obtiene la respuesta plástica perfecta, sin embargo, presenta oscilaciones en el campo de las presiones, tal como se aprecia en la figura 5f.

Figura 5. Comparación de los campos de presión y deformación plástica equivalente obtenidas con MEX-FEM y formulación estándar en desplazamientos.

Finalmente, la figura 6 compara la curva Desplazamiento-Reacción obtenida con la formulación irreducible en desplazamientos y la MEX-FEM. Nótese cómo MEX-FEM captura satisfactoriamente el comportamiento plástico perfecto y la carga última, cuyo valor predicho es de 15.16 y 15.22 para los modelos de MC y DP, respectivamente, en excelente acuerdo con la solución teórica. El efecto de bloqueo de la solución obtenida con elemento lineal de la formulación irreducible (FI − P1) en ambos modelos de plasticidad se puede apreciar por la inexistencia de una carga límite en el rango perfectamente plástico. El elemento cuadrático de la formulación irreducible (FI − P2) captura el comportamiento plástico perfecto, pero la carga última predicha es de 15.56, ligeramente por encima del valor teórico.

Figura 6. Curva Desplazamiento-Reacción.

4.2. Barra con agujero a tracción. Localización en tensión y deformación plana

Este ejemplo consiste en una barra sometida a tracción en sus extremos libres. El objetivo de este problema es doble: 1) determinar numéricamente el ángulo de localización θ_loc variando el ángulo de fricción interna ϕ y 2) investigar el efecto del coeficiente de Poisson en la dirección de la banda de localización. Para el caso de plasticidad incompresible (ϕ = 0^o), se toma c_u = 5.00. Para los demás casos c_u = 1.

Dada la simetría del problema, sólo se discretiza la parte superior derecha. El dominio computacional se discretiza en 3758 nodos y 7274 elementos. El problema se analiza tanto para estados de tensión y deformación plana. La figura 7 muestra los datos geométricos y la malla de elementos finitos empleada. Las dimensiones l, b y d se toman como 40 m, 20 m y 2 m respectivamente. El espesor se toma como 1 m. Las propiedades del material son: módulo de Young E = 10⁷ Pa, coeficiente de Poisson ν = {0, 0.15, 0.30}, cohesión inicial c₀ = 10⁴ Pa, ángulo de fricción ϕ = {0^o, 15^o, 30^o, 45^o, 60^o}.

Figura 7. Geometría y malla empleada.

Se considera primero el caso de deformación plana. Los resultados de los análisis se presentan en la tabla 3. Es notable la precisión obtenida por MEX-FEM para capturar la correcta dirección de la banda de localización en ambos modelos de MC y DP y el excelente acuerdo con la solución analítica. La figura 8 muestra el mecanismo de fallo para diferentes valores del ángulo de fricción interna ϕ en un modelo en deformación plana de MC predicho por MEX-FEM. Los resultados obtenidos son óptimos ya que ĺa localización pasa a través de una banda de elementos y está virtualmente libre de la dependencia de la malla. Se obtienen resultados similares en el modelo de DP.

Tabla 3. Ángulo de localización θ_loc para el caso de tracción uniaxial pura en deformación plana.

Ángulo de fricción ϕ	θloc analítico DP	θloc numérico DP	θloc analítico MC	θloc numérico MC
00	45o	44 . 12o	45o	44 . 54o
150	36 . 40o	36 . 35o	37 . 5o	37 . 38o
300	28 . 15o	29 . 65o	30o	30 . 25o
450	20 . 44o	23 . 29o	22 . 5o	23 . 29o
600	13 . 28o	15 . 52o	15o	15 . 52o

Figura 8. Deformación plástica equivalente para diferentes valores de ϕ en modelo MC en deformación plana.

La tabla 4 compara los ángulos de localización obtenidos para el modelo de DP en tensión plana. Nuevamente el ángulo predicho por la formulación propuesta concuerda con la solución teórica. La figura 9 muestra el campo de las deformaciones plásticas equivalentes obtenido.

Tabla 4. Ángulo de localización θ_loc para el caso de tracción uniaxial pura en tensión plana en un modelo DP

Ángulo de fricción ϕ	θloc analítico DP	θloc numérico DP
00	35 . 26o	35 . 25o
150	28 . 98o	30 . 06o
300	22 . 66o	24 . 07o
450	16 . 57o	18 . 88o

Figura 9. Deformación plástica equivalente para diferentes valores de ϕ en modelo DP.

La figura 10 muestra la banda de localización obtenida para el modelo de MC en un estado de tensión plana. Debido a la singularidad que existe en la superficie de MC cuando el estado de tensión es tracción pura, se aplica una tensión de compresión t_x < 0 y una tensión de tracción t_x > 0 en la dirección x de magnitud 1000 Pa, de tal manera que la dirección del flujo quede definida en cada uno de los planos A − C o A − B (véase fig. 1b). Para el estado en que t_x < 0 se obtiene la respuesta estándar del fallo de MC. Sin embargo, en el caso t_x > 0 se obtiene una localización horizontal de deformaciones, el mismo comportamiento que se obtendría con un modelo de plasticidad de Rankine. La tabla 5 compara los resultados obtenidos con MEX-FEM y la solución teórica.

Figura 10. Deformación plástica equivalente para diferentes valores de ϕ en modelo MC en tensión plana.

Tabla 5. Ángulo de localización θ_loc para t_x ≠ 0 en tensión plana en un modelo MC

Ángulo de fricción ϕ	Estado	θloc analítico MC	θloc numérico MC
00	ty > 0 y tx > 0	0o	0 . 00o
00	ty > 0 y tx < 0	45o	44 . 95o
300	ty > 0 y tx > 0	0o	0 . 00o
300	ty > 0 y tx < 0	30o	30 . 84o

La tabla 6 compara los resultados de la banda de localización obtenido por MEX-FEM para diferentes coeficientes de Poisson. Se aprecia que el ángulo de localización depende únicamente del vector de flujo y el estado tensional, tal y como predice el resultado analítico.

Tabla 6. Ángulo de localización θ_loc para el caso uniaxial puro t_y ≠ 0 y t_x = 0 en deformación plana para valores distintos de coeficiente de Poisson

ϕ	ν	θloc analítico DP	θloc numérico DP	θloc	θloc numérico MC
00	0.0	45o	44 . 94o	45o	44 . 94o
00	0.15	45o	44 . 12o	45o	44 . 94o
300	0.0	28 . 15o	29 . 65O	30o	29 . 65o
300	0.15	28 . 15o	29 . 65o	30o	29 . 65o

Finalmente, las figuras 11 muestran las curvas Desplazamiento-Reacción obtenidas para el modelo de MC en un estado de deformación plana y para el modelo de DP en tensión plana, respectivamente. Obsérvese que la carga última decrece a medida que el ángulo de fricción interna ϕ aumenta.

Figura 11. Curvas Desplazamiento-Reacción.

4.3. Cilindro de pequeño espesor con agujero a tracción y torsión longitudinales. Localización tridimensional

El último problema presentado es un cilindro hueco de pequeño espesor sometido a dos solicitaciones distintas: tracción y torsión según el eje longitudinal del cilindro. Las dimensiones del cilindro son: h = 1.95m, radio externo r = 0.50m y espesor t = 0.05m; el cilindro presenta un hueco rectangular en su centro a fin de provocar una concentración de tensiones y la formación de una banda helicoidal de localización de deformación plástica. Para esta geometría y las solicitaciones estudiadas el estado resultante es de tensión plana, por lo que los ángulos de localización son los de la tabla 2 ( [7]; [10] ;  [34]). Dada la simetría de la geometría, sólo se analiza la parte superior. El modelo geométrico y la discretización estructurada empleada, con 15357 nodos y 59880 elementos tetraedros MEX-FEM P1-P1, se muestran en la figura 12. En las simulaciones se usan las propiedades materiales: densidad ρ = 100 Kg/m³, módulo de Young E = 10⁵Pa, coeficiente de Poisson ν = 0.30, límite elástico f_y = c₀ = 400 Pa   y energía de fractura ${\textstyle G_{f}=5}$ N/m. Tanto la tracción como la torsión se aplican mediante desplazamientos impuestos en la superficie superior a una velocidad de 0.001 m/s.

Figura 12. Geometría y malla empleada.

En primer lugar, se estudia el caso de tracción longitudinal. Se utiliza como modelo de plasticidad el modelo de von Mises (Drucker-Prager sin fricción). En la figura 13 puede observarse como se forman dos líneas de deslizamiento en forma de espirales simétricas que se inician en el hueco rectangular y se propagan hacia los extremos del cilindro a ±45^o con el plano horizontal, con la dirección de la tensión principal menor, ([7]; [10] ;  [34]). La solución numérica coincide exactamente con la solución analítica.

Figura 13. Deformaciones plásticas equivalentes en el modelo de VM en tracción pura. Ángulo de localización teórico: ${\textstyle {\theta }_{loc}^{teo}=\pm 45.0^{o}}$ con el eje vertical. Ángulo de localización numérico: ${\textstyle {\theta }_{loc}^{num}=\pm 45.0^{o}}$.

En segundo lugar, se estudia el caso de torsión longitudinal. Se utiliza como modelo de plasticidad el modelo de Mohr-Coulomb, con un ángulo de fricción de 45^o. En este caso, el estado tensional en las paredes del cilindro es de cortante puro, con tensiones principales iguales y de signos opuestos actuando a 45^o con el plano horizontal. Según los resultados analíticos en la tabla 2 ([7]; [10] ;  [34]), las líneas de deslizamiento se forman a ±22 . 5^o con la dirección principal menor, es decir, a 22 . 5^o y 67 . 5^o con el plano horizontal. Ambas soluciones alternativas se muestran en las figuras 14 y 15, respectivamente. Para obtener una u otra, se perturba ligeramente el problema de torsión pura con una tracción/compresión longitudinal, respectivamente.

Figura 14. Deformaciones plásticas equivalentes en el modelo de MC con ϕ = 45o en torsión con perturbación de tracción. Ángulo de localización teórico: ${\textstyle {\theta }_{loc}^{teo}=22.5^{o}}$ con el plano horizontal. Ángulo de localización numérico: ${\textstyle {\theta }_{loc}^{num}=23.{62}^{o}}$.

Figura 15. Deformaciones plásticas equivalentes en el modelo de MC con ϕ = 45o en torsión pura. Ángulo de localización teórico: ${\textstyle {\theta }_{loc}^{teo}=67.5^{o}}$ con el plano horizontal. Ángulo de localización numérico: ${\textstyle {\theta }_{loc}^{num}=67.{53}^{o}}$.

Tanto el caso de tracción como el de torsión, se obtienen las soluciones analíticas con independencia de la orientación de la malla. La formulación propuesta es capaz de dar una solución satisfactoria al problema tridimensional de la localización, sin necesidad de emplear un complejo algoritmo de rastreo tridimensional.

5. Conclusión

En este trabajo se presenta la aplicación de MEX-FEM en problemas de mecánica de sólidos que involucran plasticidad y localización de deformaciones. Esta formulación elude la condición de estabilidad LBB sobre las formulaciones mixtas utilizando el método de las sub-escalas ortogonales. Los ejemplos presentados en 2D y 3D demuestran que con esta formulación se obtienen campos de deformaciones y tensiones, capturando correctamente los mecanismos de colapso, las cargas últimas y prediciendo correctamente la orientación de la banda de localización sin necesidad de emplear un algoritmo de rastreo.

Agradecimientos

Nelson Lafontaine agradece a la Agencia Española De Cooperación Internacional Para El Desarrollo (AECID) por el soporte económico dado en el programa de Becas MAEC-AECID. El autor también agradece a los profesores M.W Yuan y Chen Pu por su incondicional apoyo. El trabajo se enmarca en el The Seventh Framework Programme (FP7/2007-2013) del ERC bajo el acuerdo n^o 611636 (NUMEXAS) y en el “Excellence Programme for Knowledge Generation by MINECO” a través del proyecto EACY (Enhanced accuracy computational and experimental framework for strain localization and failure mechanisms, ref. MAT2013-48624-C2-1-P).

References

1. [1] I.B. Ronaldo; A finite element model for strain localization analysis of strongly discontinuous fields based on standard galerkin approximation; Computer Methods in Applied Mechanics and Engineering, 190 (11-12) (2000), pp. 1529–1549
2. [2] Ivo Babuska; Error-bounds for finite element method; Numerische Mathematik, 16 (4) (1971), pp. 322–333
3. [3] D. Bigoni, T. Hueckel; Uniqueness and localization-i. associative and non-associative elastoplasticity; International Journal of Solids and Structures, 28 (2) (1991), pp. 197–213
4. [4] M. Cervera, M. Chiumenti; Size effect and localization in {J2} plasticity; International Journal of Solids and Structures, 46 (17) (2009), pp. 3301–3312
5. [5] M. Cervera, M. Chiumenti, C.Agelet de Saracibar; Softening, localization and stabilization: capture of discontinuous solutions in j2 plasticity; International Journal for Numerical and Analytical Methods in Geomechanics, 28 (5) (2004), pp. 373–393
6. [6] M. Cervera, M. Chiumenti, L. Benedetti, R. Codina; Mixed stabilized finite element methods in nonlinear solid mechanics. part iii: Compressible and incompressible plasticity; Computer Methods in Applied Mechanics and Engineering, 285 (0) (2015), pp. 752–775
7. [7] M. Cervera, M. Chiumenti, D. Di Capua; Benchmarking on bifurcation and localization in {J2} plasticity for plane stress and plane strain conditions; Computer Methods in Applied Mechanics and Engineering, 241-244 (0) (2012), pp. 206–224
8. [8] M. Cervera, M. Chiumenti, R. Codina; Mixed stabilized finite element methods in nonlinear solid mechanics: Part i: Formulation; Computer Methods in Applied Mechanics and Engineering, 199 (37-40) (2010), pp. 2559–2570
9. [9] M. Cervera, M. Chiumenti, R. Codina; Mixed stabilized finite element methods in nonlinear solid mechanics: Part ii: Strain localization; Computer Methods in Applied Mechanics and Engineering, 199 (37-40) (2010), pp. 2571–2589
10. [10] M. Cervera, J.-Y. Wu; On the conformity of strong, regularized, embedded and smeared discontinuity approaches for the modeling of localized failure in solids; International Journal of Solids and Structures, 71 (2015), pp. 19–38
11. [11] R. Chambon, S. Crochepeyre, J. Desrues; Localization criteria for non-linear constitutive equations of geomaterials; Mechanics of Cohesive-frictional Materials, 5 (1) (2000), pp. 61–82
12. [12] R. Codina; Stabilization of incompresssibility and convection through orthogonal sub-scales in finite elements methods; Comput. Meth Appl Mech Eng, 190 (2000), pp. 1579–1599
13. [13] R. Codina; Stabilized finite element approximation of transient incompressible flows using orthogonal subscales; Computer Methods in Applied Mechanics and Engineering, 191 (39-40) (2002), pp. 4295–4321
14. [14] R. Codina; Analysis of a stabilized finite element approximation of the oseen equations using orthogonal subscales; Appl. Numer. Math., 58 (3) (mar 2008), pp. 264–283
15. [15] P. Dadvand, R. Rossi, M. Gil, X. Martorell, J. Cotela, E. Juanpere, S.R. Idelsohn, E. Oñate; Migration of a generic multi-physics framework to hpc environments; Computers & Fluids, 80 (0) (2013), pp. 301–309
16. [16] P. Dadvand, R. Rossi, E. Oñate; An object-oriented environment for developing finite element codes for multi-disciplinary applications; Archives of Computational Methods in Engineering, 17 (3) (2010), pp. 253–297
17. [17] C. Agelet de Saracibar, M. Chiumenti, Q. Valverde, M. Cervera; On the orthogonal subgrid scale pressure stabilization of finite deformation {J2} plasticity; Computer Methods in Applied Mechanics and Engineering, 195 (9-12) (2006), pp. 1224–1251
18. [18] F. Brezzi, M. Fortin, D. Marini; Mixed finite element methods; Springer (1991)
19. [19] GiD. The personal pre and post processor. http://gid.cimne.upc.es. 2009.
20. [20] T.J.R. Hughes, G. Scovazzi, L.P. Franca; Multiscale and Stabilized Methods; John Wiley & Sons, Ltd (2004)
21. [21] Thomas J.R. Hughes; Multiscale phenomena: Greens functions, the dirichlet-to-neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods; Computer Methods in Applied Mechanics and Engineering, 127 (1-4) (1995), pp. 387–401
22. [22] M. Iordache, K. Willam; Localized failure analysis in elastoplastic cosserat continua; Computer Methods in Applied Mechanics and Engineering, 151 (3-4) (1998), pp. 559–586 Containing papers presented at the Symposium on Advances in Computational Mechanics
23. [23] N.M. Lafontaine, R. Rossi, M. Cervera, M. Chiumenti; Explicit mixed strain-displacement finite element for dynamic geometrically non-linear solid mechanics; Computational Mechanics (2015), pp. 1–17
24. [24] N.M. Lafontaine, R. Rossi, M. Cervera, M. Chiumenti; Formulación mixta estabilizada explícita de elementos finitos para sólidos compresibles y quasi-incompresibles; Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingenier’í (2015)
25. [25] Y. Leroy, M. Ortiz; Finite element analysis of transient strain localization phenomena in frictional solids; International Journal for Numerical and Analytical Methods in Geomechanics, 14 (2) (1990), pp. 93–124
26. [26] E.d.S. Neto, D.R.J. Peric; Computational Methods for Plasticity; Wiley (2008)
27. [27] N. Saabye Ottosen, K. Runesson; Properties of discontinuous bifurcation solutions in elasto-plasticity; International Journal of Solids and Structures, 27 (4) (1991), pp. 401–421
28. [28] E. Oñate, J. Rojek, R.L. Taylor, O.C. Zienkiewicz; Finite calculus formulation for incompressible solids using linear triangles and tetrahedra; International Journal for Numerical Methods in Engineering, 59 (11) (2004), pp. 1473–1500
29. [29] J.W. Rudnicki, J.R. Rice; Conditions for the localization of deformation in pressure-sensitive dilatant materials; Journal of the Mechanics and Physics of Solids, 23 (6) (1975), pp. 371–394
30. [30] K. Runesson, N. Saabye Ottosen, P. Dunja; Discontinuous bifurcations of elastic-plastic solutions at plane stress and plane strain; International Journal of Plasticity, 7 (1-2) (1991), pp. 99–121
31. [31] J.C. Simo, T.J.R. Hughes; Computational Inelasticity; Springer-Verlag, New-York (1998)
32. [32] P. Steinmann, K. Willam; Finite element analysis of elastoplastic discontinuities; Journal of Engineering Mechanics, 120 (11) (1994), pp. 2428–2442
33. [33] Q. Valverde, M. Chiumenti, M. Cervera, C.Agelet de Saracibar; Formulación estabilizada de elementos finitos triangulares y tetraédricos para problemas de incompresibilidad en pequeñas deformaciones; Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingenier’í, 21 (4) (2005)
34. [34] J.-Y. Wu, M. Cervera; On the equivalence between traction- and stress-based approaches for the modeling of localized failure in solids; Journal of the Mechanics and Physics of Solids, 82 (2015), pp. 137–163