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 [[#bib0115|[23]]] ;  [[#bib0120|[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)[[#bib0100|[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 [[#bib0030|[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 [[#bib0020|[4]]]; [[#bib0025|[5]]]; [[#bib0035|[7]]]; [[#bib0085|[17]]] ;  [[#bib0165|[33]]] han empleado esta formulación para la solución de problemas de plasticidad ''J''2 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 [[#bib0040|[8]]] ;  [[#bib0045|[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 [[#bib0030|[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 [[#sec0010|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 [[#sec0030|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 [[#sec0050|4]] presenta una serie de ejemplos numéricos que validan el buen comportamiento de los elementos de la formulación propuesta.

<span id='sec0010'></span>
==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 <math display="inline">\overline{\boldsymbol{t}}=\boldsymbol{\sigma }\boldsymbol{n}:\partial {\Omega }_{\sigma }\rightarrow {\mathbb{R}}^{n_{dim}}</math>, los desplazamientos '''''u''''' y aceleraciones <math display="inline">\ddot{\boldsymbol{u}}</math> en ∂Ω<sub>''u''</sub> y las fuerzas por unidad de volumen <math display="inline">\boldsymbol{b}:\Omega \rightarrow {\mathbb{R}}^{n_{dim}}</math>, hallar el campo de los desplazamientos, velocidades, aceleraciones y deformaciones en cualquier instante de tiempo <math display="inline">t\in I</math> siendo <math display="inline">I=(0,T)</math> el intervalo de tiempo de interés, tales que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\nabla \cdot \boldsymbol{\sigma }+\boldsymbol{b}=\rho \ddot{\boldsymbol{u}}\quad en\quad \Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

donde <math display="inline">\Omega \subset {\mathbb{R}}^{n_{dim}}</math> es el volumen ocupado por el sólido en un espacio de ''n''<sub>''dim''</sub> 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 '''''ɛ'''''<sub>''e''</sub> y otra plástica '''''ɛ'''''<sub>''p''</sub>,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\varepsilon }={\boldsymbol{\varepsilon }}_e+</math><math>{\boldsymbol{\varepsilon }}_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

La ecuación constitutiva en plasticidad es:

<span id='eq0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\sigma }=\boldsymbol{C}_0:\left(\boldsymbol{\varepsilon }-\right. </math><math>\left. {\boldsymbol{\varepsilon }}_p\right)=\boldsymbol{C}_0:{\boldsymbol{\varepsilon }}_e</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

donde '''''C'''''<sub>0</sub> es el tensor constitutivo elástico inicial. El problema queda definido formulando una ley de evolución del tipo <math display="inline">\dot{{\boldsymbol{\varepsilon }}_p}={\dot{\boldsymbol{\varepsilon }}}_p(\boldsymbol{\sigma })</math> para las deformaciones plásticas.

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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 }=</math><math>\boldsymbol{C}(\boldsymbol{\varepsilon },{\boldsymbol{\varepsilon }}_p):\boldsymbol{\varepsilon }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

siendo '''''C'''''('''''ɛ''''', '''''ɛ'''''<sub>''p''</sub>) 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 [[#bib0030|[6]]] y [[#bib0115|[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 [[#eq0010|(1)]] y [[#eq0015|(2)]], las variables del problema '''''u''''' y '''''ɛ''''' están sujetas a unas condiciones iniciales en ''t'' = ''t''<sub>0</sub>, es decir; <math display="inline">\boldsymbol{u}\vert _{t=t_0}=\boldsymbol{u}_0</math>, <math display="inline">\dot{\boldsymbol{u}}\vert _{t=t_0}=\boldsymbol{v}_0</math> y <math display="inline">\boldsymbol{\varepsilon }\vert _{t=t_0}={\boldsymbol{\varepsilon }}_0=</math><math>{\nabla }^s\boldsymbol{u}_0</math>. Después de multiplicar las ecuaciones [[#eq0010|(1)]] y [[#eq0015|(2)]] por sus respectivas funciones de peso e integrando por partes la ecuación [[#eq0010|(1)]], la formulación variacional del problema es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\nabla }^s\boldsymbol{\omega }:\boldsymbol{\sigma }\quad d\Omega +</math><math>{\int }_{\Omega }\boldsymbol{\omega }\cdot \rho \ddot{\boldsymbol{u}}\quad d\Omega =</math><math>{\int }_{\partial \Omega }\boldsymbol{\omega }\cdot \overline{\boldsymbol{t}}\quad d\Gamma +</math><math>{\int }_{\Omega }\boldsymbol{\omega }\cdot \boldsymbol{b}\quad d\Omega \quad \quad \forall \boldsymbol{\omega }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

donde <math display="inline">\boldsymbol{\omega }\in V</math> y <math display="inline">\boldsymbol{\gamma }\in T</math> son las funciones de prueba para el campo de los desplazamientos y el campo de las deformaciones, respectivamente; <math display="inline">V</math> y <math display="inline">T</math> son los espacios de los desplazamientos y deformaciones admisibles.

Para la definición discreta del problema, se introducen en las ecuaciones [[#eq0040|(6)]] y [[#eq0045|(7)]] un campo de desplazamientos discretos, '''''u''''' = '''''u'''''<sub>''h''</sub>, aceleraciones discretas <math display="inline">\ddot{\boldsymbol{u}}={\ddot{\boldsymbol{u}}}_h</math> y deformaciones discretas '''''ɛ''''' = '''''ɛ'''''<sub>''h''</sub> tales que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:\boldsymbol{\sigma }\quad d\Omega +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho {\ddot{\boldsymbol{u}}}_h\quad d\Omega =</math><math>{\int }_{\partial \Omega }{\boldsymbol{\omega }}_h\cdot \overline{\boldsymbol{t}}\quad d\Gamma +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \boldsymbol{b}\quad d\Omega \quad \quad \forall {\boldsymbol{\omega }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

en donde <math display="inline">{\boldsymbol{\omega }}_h\in V_h</math> y <math display="inline">{\boldsymbol{\gamma }}_h\in T_h</math> son las funciones de interpolación definidas sobre el espacio de los elementos finitos <math display="inline">V_h</math> y <math display="inline">T_h</math>, 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 [[#bib0010|[2]]] ;  [[#bib0090|[18]]], por lo que la solución obtenida con las ecuaciones [[#eq0055|(8)]] y [[#eq0060|(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)[[#bib0105|[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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{u}=\boldsymbol{u}_h+\tilde{\boldsymbol{u}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

Asimismo, derivando en el tiempo,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\dot{\boldsymbol{u}}={\dot{\boldsymbol{u}}}_h+\dot{\tilde{\boldsymbol{u}}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

donde <math display="inline">{\boldsymbol{\varepsilon }}_h\in T_h</math> y <math display="inline">(\boldsymbol{u}_h,{\dot{\boldsymbol{u}}}_h,{\ddot{\boldsymbol{u}}}_h)\in V_h</math> son las deformaciones, desplazamientos, velocidades y aceleraciones en la escala gruesa; <math display="inline">\tilde{\boldsymbol{\varepsilon }}\in \tilde{T}</math> y <math display="inline">(\tilde{\boldsymbol{u}},\dot{\tilde{\boldsymbol{u}}},\ddot{\tilde{\boldsymbol{u}}})\in \tilde{V}</math> son las deformaciones y desplazamientos (con sus derivadas temporales) de la sub-escala y <math display="inline">(T_h,V_h)</math> e <math display="inline">(\tilde{T},\tilde{V})</math> 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:

<span id='eq0095'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\sigma }=\boldsymbol{C}_0:\left({\boldsymbol{\varepsilon }}_h+\right. </math><math>\left. \tilde{\boldsymbol{\varepsilon }}\right)=\boldsymbol{C}_0:{\boldsymbol{\varepsilon }}_h+</math><math>\boldsymbol{C}_0:\tilde{\boldsymbol{\varepsilon }}=</math><math>{\boldsymbol{\sigma }}_h+\tilde{\boldsymbol{\sigma }}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Introduciendo las ecuaciones [[#eq0070|(10)]], [[#eq0085|(12)]], [[#eq0090|(13)]], [[#eq0075|(11)]] y [[#eq0095|(14)]] en las ecuaciones [[#eq0040|(6)]] y [[#eq0045|(7)]] se obtiene:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:{\boldsymbol{\sigma }}_h\quad d\Omega +</math><math>{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:\tilde{\boldsymbol{\sigma }}\quad d\Omega +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho \left({\ddot{\boldsymbol{u}}}_h+\right. </math><math>\left. \ddot{\tilde{\boldsymbol{u}}}\right)\quad d\Omega =</math><math>F_{ext}({\boldsymbol{\omega }}_h)\quad \forall {\boldsymbol{\omega }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

donde <math display="inline">F_{ext}({\boldsymbol{\omega }}_h)={\int }_{\partial \Omega }{\boldsymbol{\omega }}_h\cdot \overline{\boldsymbol{t}}\quad d\Gamma +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \boldsymbol{b}\quad d\Omega </math>. Si bien las ecuaciones [[#eq0105|(15)]] y [[#eq0115|(17)]] se plantean y se resuelven en el espacio de los elementos finitos, no así las ecuaciones [[#eq0110|(16)]] y [[#eq0120|(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 [[#eq0105|(15)]] y [[#eq0115|(17)]]. Aquí se emplea el método de las Sub-escalas Ortogonales OSS.

Las ecuaciones [[#eq0105|(15)]] y [[#eq0110|(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 [[#bib0120|[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 <math display="inline">\parallel \tilde{\boldsymbol{\varepsilon }}\parallel <<\parallel {\boldsymbol{\varepsilon }}_h\parallel </math> (véase [[#bib0030|[6]]]; [[#bib0040|[8]]]; [[#bib0115|[23]]] ;  [[#bib0120|[24]]]), es posible aproximar

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\varepsilon }}_p={\boldsymbol{\varepsilon }}_p(\boldsymbol{\sigma })\simeq {\boldsymbol{\varepsilon }}_p({\boldsymbol{\sigma }}_h)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

lo que implica que las tensiones continuas son:

<span id='eq0130'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\sigma }=\boldsymbol{C}_0:\left({\boldsymbol{\varepsilon }}_h+\right. </math><math>\left. \tilde{\boldsymbol{\varepsilon }}-{\boldsymbol{\varepsilon }}_p({\boldsymbol{\sigma }}_h)\right)=</math><math>{\boldsymbol{\sigma }}_h+\tilde{\boldsymbol{\sigma }}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 [[#eq0130|(20)]] es análoga a la expresión [[#eq0095|(14)]]. De esta manera se tiene el mismo formalismo y procedimiento que en el caso elástico [[#bib0120|[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) [[#bib0060|[12]]]; [[#bib0065|[13]]] ;  [[#bib0070|[14]]]. Por practicidad se introduce el operador de proyección ortogonal <math display="inline">P_h^{\perp }(\bullet )=(\bullet )-P_h(\bullet )</math>. 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 <math display="inline">\tilde{\xi }\in [0,1]</math> se tiene que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tilde{\boldsymbol{u}}}^{n+1}\simeq {\tau }_{u_t}\frac{\rho }{\Delta t^2}\left((2-\right. </math><math>\left. \tilde{\xi }){\tilde{\boldsymbol{u}}}^n-(1-\right. </math><math>\left. \tilde{\xi }){\tilde{\boldsymbol{u}}}^{n-1}\right)+</math><math>{\tau }_{u_t}P^{\perp }\left(\nabla \cdot {\boldsymbol{\sigma }}_h^n\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

Los términos <math display="inline">{\tau }_{u_t}</math> y ''τ''<sub>''ɛ''</sub> son los parámetros de estabilización dados por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tau }_{u_t}={\left(\frac{\rho }{\Delta t^2}+B(\mu )\frac{{\mu }_0}{c_uh_e^2}\right)}^{-1}\quad y\quad {\tau }_{\varepsilon }=</math><math>c_{\varepsilon }B(\mu )\frac{h_e}{L_0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

donde ''c''<sub>''u''</sub> > 0 y ''c''<sub>''ɛ''</sub> > 0 son constantes algorítmicas adimensionales, ''h''<sub>''e''</sub> es la longitud del elemento finito, ''L''<sub>0</sub> es la longitud característica del problema, ''μ''<sub>0</sub> = 2''G'' es el módulo inicial de rigidez al corte y ''β''(''μ'') = ''μ''/''μ''<sub>0</sub>, siendo <math display="inline">\mu =\frac{\parallel \boldsymbol{S}_h\parallel }{\parallel dev({\boldsymbol{\varepsilon }}_h)\parallel }</math> el módulo de corte secante [[#bib0140|[28]]] ;  [[#bib0165|[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 <math display="inline">\overline{\mu }=(1-\overline{\eta }){\mu }_{old}+\overline{\eta }{\mu }_{new}</math> donde <math display="inline">\overline{\eta }\in \lbrace 0,1\rbrace </math> 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 <math display="inline">p=\frac{1}{3}tr(\boldsymbol{\sigma })</math> en situaciones cuasi-incompresibles. El término que estabiliza el campo de la tensión media es el gradiente ∇''p''<sub>''h''</sub> contenido dentro de ∇ · '''''σ'''''<sub>''h''</sub> = ∇ · '''''S'''''<sub>''h''</sub> + ∇ ''p''<sub>''h''</sub> siendo '''''S'''''<sub>''h''</sub> el tensor desviador de tensiones. Se ha demostrado efectivo sustituir el último término de la ecuación [[#eq0140|(21)]], <math display="inline">P_h^{\perp }(\nabla \cdot {\boldsymbol{\sigma }}_h)</math>, por <math display="inline">P_h^{\perp }(\nabla p_h)</math>, quedando

<span id='eq0155'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tilde{\boldsymbol{u}}}^{n+1}\simeq {\tau }_{u_t}\frac{\rho }{\Delta t^2}\left((2-\right. </math><math>\left. \tilde{\xi }){\tilde{\boldsymbol{u}}}^n-(1-\right. </math><math>\left. \tilde{\xi }){\tilde{\boldsymbol{u}}}^{n-1}\right)+</math><math>{\tau }_{u_t}\left(\nabla p_h^n-P_h\left(\nabla p_h^n\right)\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

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

Las funciones de prueba '''''ω'''''<sub>''h''</sub> para el espacio de los desplazamientos y '''''γ'''''<sub>''h''</sub> para el espacio de las deformaciones, en las ecuaciones [[#eq0105|(15)]] y [[#eq0115|(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 [[#bib0120|[24]]]:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{u}_h=\sum_{i=1}^4\boldsymbol{N}_u^i\boldsymbol{u}_h^i\quad y\quad {\boldsymbol{\varepsilon }}_h=</math><math>\sum_{i=1}^4\boldsymbol{N}_{\boldsymbol{\varepsilon }}^i{\boldsymbol{\varepsilon }}_h^i</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

siendo <math display="inline">\boldsymbol{N}_u^i</math> la submatriz diagonal de 3 × 3 <math display="inline">(\boldsymbol{N}_u^i=diag\lbrace N_i,N_i,N_i\rbrace )</math> de las funciones de interpolación del campo de los desplazamientos en el nodo ''i  ''<math display="inline">(\boldsymbol{u}_h^i=\lbrace u_{x,h}^i,u_{y,h}^i,u_{z,h}^i\rbrace )</math> y <math display="inline">\boldsymbol{N}_{\boldsymbol{\varepsilon }}^i</math> la submatriz diagonal de 6 × 6 <math display="inline">(\boldsymbol{N}_{\boldsymbol{\varepsilon }}^i=diag\lbrace N_i,N_i,N_i,N_i,N_i,N_i\rbrace )</math> de las funciones de interpolación del campo de las deformaciones en el nodo ''i  ''<math display="inline">({\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 )</math>. Por otro lado, en la ecuación [[#eq0105|(15)]] aparece el término <math display="inline">{\nabla }^s{\boldsymbol{\omega }}_h=\boldsymbol{B}_u^T</math> que no es más que el clásico operador matricial de gradiente simétrico discreto, en el cual '''''B'''''<sub>''u''</sub> = ['''''B'''''<sub>1</sub>'''''B'''''<sub>2</sub>'''''B'''''<sub>3</sub>'''''B'''''<sub>4</sub>], donde cada sub-matriz '''''B'''''<sub>''i''</sub> se escribe como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

Asimismo, el operador de divergencia ∇ · '''''γ'''''<sub>''h''</sub> en la ecuación [[#eq0115|(17)]] toma la forma de <math display="inline">\boldsymbol{B}_{\varepsilon }=[\boldsymbol{B}_1^T\boldsymbol{B}_2^T\boldsymbol{B}_3^T\boldsymbol{B}_4^T]</math>.

El término inercial de la ecuación [[#eq0105|(15)]] da lugar al vector de fuerzas de inercia, <math display="inline">\boldsymbol{M}{\ddot{\boldsymbol{u}}}_h</math>, donde <math display="inline">\boldsymbol{M}={\int }_{{\Omega }_e}\boldsymbol{N}_u^T\rho \boldsymbol{N}_ud\Omega </math> 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 [[#bib0115|[23]]] ;  [[#bib0120|[24]]]:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\varepsilon }}_h^n=P_h\left({\nabla }^s\boldsymbol{u}^n\right)-</math><math>P_h\left(\left(\nabla \cdot {\boldsymbol{\gamma }}_h\right)\cdot {\tilde{\boldsymbol{u}}}^n\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\ddot{\boldsymbol{u}}}_h^n=\frac{1}{\Delta t^2}\left(\boldsymbol{u}_h^{n+1}-\right. </math><math>\left. 2\boldsymbol{u}_h^n+\boldsymbol{u}_h^{n-1}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

El algortimo de resolución de las ecuaciones mixta explícita sigue los pasos siguientes [[#bib0120|[24]]]:
* 1. Cómputo de las fuerzas internas estabilizadas <math display="inline">\boldsymbol{F}_{int}^{stab,n}</math>.
* 2. Cómputo de los desplazamientos en ''t'' = ''t''<sup>''n''+1</sup>:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{u}}^{n+1}={\left[2\boldsymbol{M}+\Delta t\boldsymbol{D}\right]}^{-1}[(4\boldsymbol{M}\boldsymbol{\mbox{u}}^n-</math><math>[2\boldsymbol{M}-\Delta t\boldsymbol{D}]\boldsymbol{\mbox{u}}^{n-1}+</math><math>2\Delta t^2(\boldsymbol{F}_{ext}^n-\boldsymbol{F}_{int}^{stab,n})]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}
* 3. Evaluación de las sub-escalas de los desplazamientos <math display="inline">{\tilde{\boldsymbol{u}}}^{n+1}</math> empleando la ecuación [[#eq0155|(24)]].
* 4. Actualización de las deformaciones <math display="inline">{\boldsymbol{\varepsilon }}_h^{n+1}</math> con la ecuación [[#eq0175|(27)]] en ''t''<sup>''n''+1</sup>.
* 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 [[#bib0115|[23]]].

<span id='sec0030'></span>
==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 ''σ''<sub>''n''</sub> alcanzan un valor crítico:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\vert \tau \vert =c-{\sigma }_ntan\phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

siendo ''c'' = ''c''('''''ɛ'''''<sub>''p''</sub>) ≥ 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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Phi \left(\tau ,{\sigma }_n,c({\boldsymbol{\varepsilon }}_p),\phi \right)=</math><math>\vert \tau \vert -c({\boldsymbol{\varepsilon }}_p)+</math><math>{\sigma }_ntan\phi =0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\Phi }_1\left({\sigma }_1,{\sigma }_3,c({\boldsymbol{\varepsilon }}_p),\phi \right)=</math><math>({\sigma }_1-{\sigma }_3)+({\sigma }_1+{\sigma }_3)sin\phi -</math><math>2c({\boldsymbol{\varepsilon }}_p)cos\phi =0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

siendo ''σ''<sub>1</sub> ≥ ''σ''<sub>2</sub> ≥ ''σ''<sub>3</sub> las tensiones principales. En el espacio de HW, el modelo de MC es una pirámide hexagonal cuyo ápice se localiza en <math display="inline">\sqrt{3}ccot\phi </math> tal como se muestra en la [[#fig0005|figura 1]].

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

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr1.jpg|center|471px|Modelo de Mohr-Coulomb.]]


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

Figura 1.

Modelo de Mohr-Coulomb.

</span>
|}

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 <math display="inline">J_2=\frac{1}{2}\boldsymbol{S}:\boldsymbol{S}</math> y la tensión hidrostática alcanzan una combinación crítica:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Phi \left(\boldsymbol{\sigma },c({\boldsymbol{\varepsilon }}_p),\phi \right)=</math><math>\sqrt{J_2(\boldsymbol{\sigma })}+\eta p-\varsigma c({\boldsymbol{\varepsilon }}_p)=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

El criterio de DP, mostrado en la [[#fig0010|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 [[#bib0130|[26]]]:

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

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr2.jpg|center|522px|Modelo de Drucker-Prager.]]


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

Figura 2.

Modelo de Drucker-Prager.

</span>
|}

<span id='eq0230'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\eta =\frac{3tan\phi }{\sqrt{9+12{tan}^2\phi }}\quad y\quad \varsigma =</math><math>\frac{3}{\sqrt{9+12{tan}^2\phi }}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 [[#bib0155|[31]]]:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\dot{\boldsymbol{\varepsilon }}}_p=\sum_{\alpha =1}^{n_{act}}{\dot{\lambda }}_{\alpha }\frac{\partial {\Phi }_{\alpha }}{\partial \boldsymbol{\sigma }}=</math><math>\sum_{\alpha =1}^{n_{act}}{\dot{\lambda }}_{\alpha }\boldsymbol{m}_{\alpha }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

donde <math display="inline">{\dot{\lambda }}_{\alpha }</math> son los parámetros de consistencia plástica y ''n''<sub>''act''</sub> son las superficies <math display="inline">J</math>''activas''. El problema consiste en determinar '''''ɛ'''''<sub>''p''</sub>, ''λ''<sub>''α''</sub> y las superficies involucradas <math display="inline">J:={\Phi }_{\alpha }\forall \alpha \in \lbrace 1,2..n_{act}\rbrace </math>. Asimismo, los parámetros de consistencia plástica siguen las condiciones complementarias de Kuhn-Tucker extendidas a modelos de plasticidad con superficies múltiples:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\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)=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

Los procedimientos de la integración constitutiva del modelo de MC y DP empleados en este trabajo están descritos en la referencia [[#bib0130|[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 <math display="inline">G_f</math> por unidad de área. En un proceso uniaxial de carga, con ablandamiento y relajación total final de las tensiones, el trabajo plástico <math display="inline">W_p</math> 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''<sub>''e''</sub> como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>W_p=\frac{G_f}{h_e}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

Suponiendo que la evolución de la cohesión sigue lo descrito en la [[#tbl0005|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''<sub>0</sub> y <math display="inline">{\overline{\boldsymbol{\varepsilon }}}_p=0</math>) hasta el desarrollo de deformaciones plásticas localizadas (''t'' =∞, ''c'' = 0 y <math display="inline">{\overline{\boldsymbol{\varepsilon }}}_p\not =0</math>) es igual a:

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Curva de ablandamiento de la cohesión

|-

! <math display="inline">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 ablandamiento\mbox{\;  }lineal\\
0 & \quad si\quad {\overline{\boldsymbol{\varepsilon }}}_p>\frac{c_0}{H_s}.
\end{array}</math>
|-

! <math display="inline">c=c_0exp\left(\frac{-2H_s}{c_0}{\overline{\boldsymbol{\varepsilon }}}_p\right)\quad para\quad 0<{\overline{\boldsymbol{\varepsilon }}}_p<\infty \quad ablandamiento\mbox{\;  }exponencial</math>
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>W_p={\int }_{t=0}^{t=\infty }{\dot{W}}_pdt={\int }_{t=0}^{t=\infty }\sigma :{\dot{\boldsymbol{\varepsilon }}}_pdt=</math><math>\frac{c_0^2}{2H_s}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

El escalar <math display="inline">H_s</math> 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 [[#bib0005|[1]]]; [[#bib0015|[3]]]; [[#bib0055|[11]]]; [[#bib0110|[22]]]; [[#bib0125|[25]]]; [[#bib0135|[27]]]; [[#bib0145|[29]]]; [[#bib0150|[30]]] ;  [[#bib0160|[32]]] han encontrado soluciones analíticas y geométricas para la orientación de las bandas de discontinuidad <math display="inline">S</math> (véase [[#fig0015|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 [[#bib0035|[7]]]; [[#bib0050|[10]]] ;  [[#bib0170|[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 [[#tbl0010|tabla 2]] muestra los valores analíticos del ángulo de localización ''θ''<sub>''loc''</sub> para los modelos de DP y MC en tensión y deformación plana. Cervera et al [[#bib0030|[6]]] han verificado numéricamente por EF los resultados analíticos de localización en un modelo de plasticidad ''J''2 y DP empleando una formulación mixta implícita.

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

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr3.jpg|center|263px|Ángulo de localización.]]


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

Figura 3.

Ángulo de localización.

</span>
|}

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

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
| <math display="inline">{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 }</math>
| <math display="inline">{tan}^2{\theta }_{loc}=\frac{{\sigma }_1-{\sigma }_2-2J_2^{1/2}\eta }{{\sigma }_1-{\sigma }_2+2J_2^{1/2}\eta }</math>
|-

!  MC
| <math display="inline">{\theta }_{loc}=\begin{array}{ll}
0^0 & si\mbox{\;  }flujo\mbox{\;  }normal\mbox{\;  }al\mbox{\;  }plano\quad A-B.\\
\pm \left({45}^o-\frac{\phi }{2}\right) & si\mbox{\;  }flujo\mbox{\;  }normal\mbox{\;  }al\mbox{\;  }plano\quad A-C.
\end{array}</math>
| <math display="inline">{\theta }_{loc}=\pm \left({45}^o-\frac{\phi }{2}\right)</math>
|}

La formulación mixta considera la ecuación geométrica [[#eq0045|(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.

<span id='sec0050'></span>
==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 <math display="inline">h_e={\left(4/\pi \cdot A_e\right)}^{\frac{1}{2}}</math> siendo ''A''<sub>''e''</sub> el área del elemento finito correspondiente. En 3''D  '' se toma <math display="inline">h_e={\left(6/\pi \cdot V_e\right)}^{\frac{1}{3}}</math> donde ''V''<sub>''e''</sub> 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''<sub>''u''</sub> = [0.25, 5] y ''c''<sub>''ɛ''</sub> = 1.0. Los análisis se han realizado con el programa de elementos finitos KRATOS [[#bib0075|[15]]] ;  [[#bib0080|[16]]], desarrollado en el Centro Internacional de Métodos Numéricos (CIMNE). Como pre y post-procesador se ha utilizado GiD [[#bib0095|[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 [[#fig0020|4]] muestran la geometría y la malla empleada en el análisis. La propiedad del material son: densidad ''ρ'' = 10<sup>4</sup> ''Kg''/''m''<sup>3</sup>, módulo de Young ''E'' = 10<sup>7</sup> ''KPa'', coeficiente de Poisson ''ν'' = 0.48, cohesión inicial ''c''<sub>0</sub> = 490 ''KPa'' y ángulo de fricción ''ϕ'' = 20<sup>''o''</sup>. Este ejemplo se analiza empleando el modelo de plasticidad perfecta de MC y DP con ''η'' y ''ς'' dados en la ecuación [[#eq0230|(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''<sub>0</sub> se toman como 5''m'' y 1''m'', respectivamente 5''m''. El dominio se discretiza empleando una malla no estructurada de 2438 nodos y 4731 elementos lineales, tanto para la formulación irreducible (''FI'' − ''P''1), como para la presente formulación (''MEX'' − ''FEMP''1''P''1). Se impone una velocidad instantánea y constante en el tiempo de 10<sup>−3</sup>''m''/''s''. Asimismo, el problema se analiza estáticamente empleando los elementos de orden cuadrático de la formulación irreducible (''FI'' − ''P''2). La solución analítica es ''P''<sub>''lim''</sub>/''c''<sub>0</sub> = 14.8, donde ''P''<sub>''lim''</sub> es la fuerza total (reacción) aplicada.

<span id='fig0020'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr4.jpg|center|391px|Geometría y malla empleada en el ensayo de Prantl.]]


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

Figura 4.

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

</span>
|}

Las figuras [[#fig0025|5]]a y [[#fig0025|5]]b 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 [[#fig0025|5]]c y [[#fig0025|5]]d muestran las distribuciones obtenidas con la formulación irreducible (''FI'' − ''P''1). 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'' − ''P''2) obtiene la respuesta plástica perfecta, sin embargo, presenta oscilaciones en el campo de las presiones, tal como se aprecia en la [[#fig0025|figura 5]]f.

<span id='fig0025'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr5.jpg|center|497px|Comparación de los campos de presión y deformación plástica equivalente ...]]


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

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.

</span>
|}

Finalmente, la [[#fig0030|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'' − ''P''1) 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'' − ''P''2) captura el comportamiento plástico perfecto, pero la carga última predicha es de 15.56, ligeramente por encima del valor teórico.

<span id='fig0030'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr6.jpg|center|349px|Curva Desplazamiento-Reacción.]]


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

Figura 6.

Curva Desplazamiento-Reacción.

</span>
|}

===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 ''θ''<sub>''loc''</sub> 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<sup>''o''</sup>), se toma ''c''<sub>''u''</sub> = 5.00. Para los demás casos ''c''<sub>''u''</sub> = 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 [[#fig0035|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<sup>7</sup> ''Pa'', coeficiente de Poisson ''ν'' = {0, 0.15, 0.30}, cohesión inicial ''c''<sub>0</sub> = 10<sup>4</sup> ''Pa'', ángulo de fricción ''ϕ'' = {0<sup>''o''</sup>, 15<sup>''o''</sup>, 30<sup>''o''</sup>, 45<sup>''o''</sup>, 60<sup>''o''</sup>}.

<span id='fig0035'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr7.jpg|center|496px|Geometría y malla empleada.]]


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

Figura 7.

Geometría y malla empleada.

</span>
|}

Se considera primero el caso de deformación plana. Los resultados de los análisis se presentan en la [[#tbl0015|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 [[#fig0040|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.

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3.

Ángulo de localización ''θ''<sub>''loc''</sub> para el caso de tracción uniaxial pura en deformación plana.

|-

! Ángulo de fricción ''ϕ''
! ''θ''<sub>''loc''</sub> analítico ''DP''
! ''θ''<sub>''loc''</sub> numérico ''DP''
! ''θ''<sub>''loc''</sub> analítico ''MC''
! ''θ''<sub>''loc''</sub> numérico ''MC''
|-

!  0<sup>0</sup>
|  45<sup>''o''</sup>
|  44 . 12<sup>''o''</sup>
|  45<sup>''o''</sup>
|  44 . 54<sup>''o''</sup>
|-

!  15<sup>0</sup>
|  36 . 40<sup>''o''</sup>
|  36 . 35<sup>''o''</sup>
|  37 . 5<sup>''o''</sup>
|  37 . 38<sup>''o''</sup>
|-

!  30<sup>0</sup>
|  28 . 15<sup>''o''</sup>
|  29 . 65<sup>''o''</sup>
|  30<sup>''o''</sup>
|  30 . 25<sup>''o''</sup>
|-

!  45<sup>0</sup>
|  20 . 44<sup>''o''</sup>
|  23 . 29<sup>''o''</sup>
|  22 . 5<sup>''o''</sup>
|  23 . 29<sup>''o''</sup>
|-

!  60<sup>0</sup>
|  13 . 28<sup>''o''</sup>
|  15 . 52<sup>''o''</sup>
|  15<sup>''o''</sup>
|  15 . 52<sup>''o''</sup>
|}

<span id='fig0040'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr8.jpg|center|px|Deformación plástica equivalente para diferentes valores de ϕ en modelo MC en ...]]


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

Figura 8.

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

</span>
|}

La [[#tbl0020|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 [[#fig0045|figura 9]] muestra el campo de las deformaciones plásticas equivalentes obtenido.

<span id='tbl0020'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 4.

Ángulo de localización ''θ''<sub>''loc''</sub> para el caso de tracción uniaxial pura en tensión plana en un modelo DP

|-

! Ángulo de fricción ''ϕ''
! ''θ''<sub>''loc''</sub> analítico ''DP''
! ''θ''<sub>''loc''</sub> numérico ''DP''
|-

!  0<sup>0</sup>
|  35 . 26<sup>''o''</sup>
|  35 . 25<sup>''o''</sup>
|-

!  15<sup>0</sup>
|  28 . 98<sup>''o''</sup>
|  30 . 06<sup>''o''</sup>
|-

!  30<sup>0</sup>
|  22 . 66<sup>''o''</sup>
|  24 . 07<sup>''o''</sup>
|-

!  45<sup>0</sup>
|  16 . 57<sup>''o''</sup>
|  18 . 88<sup>''o''</sup>
|}

<span id='fig0045'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr9.jpg|center|549px|Deformación plástica equivalente para diferentes valores de ϕ en modelo DP.]]


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

Figura 9.

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

</span>
|}

La [[#fig0050|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'''''<sub>''x''</sub> < 0 y una tensión de tracción '''''t'''''<sub>''x''</sub> > 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 [[#fig0005|fig. 1]]b). Para el estado en que '''''t'''''<sub>''x''</sub> < 0 se obtiene la respuesta estándar del fallo de MC. Sin embargo, en el caso '''''t'''''<sub>''x''</sub> > 0 se obtiene una localización horizontal de deformaciones, el mismo comportamiento que se obtendría con un modelo de plasticidad de Rankine. La [[#tbl0025|tabla 5]] compara los resultados obtenidos con MEX-FEM y la solución teórica.

<span id='fig0050'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr10.jpg|center|527px|Deformación plástica equivalente para diferentes valores de ϕ en modelo MC en ...]]


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

Figura 10.

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

</span>
|}

<span id='tbl0025'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 5.

Ángulo de localización ''θ''<sub>''loc''</sub> para '''''t'''''<sub>''x''</sub> ≠ 0 en tensión plana en un modelo MC

|-

! Ángulo de fricción ''ϕ''
! Estado
! ''θ''<sub>''loc''</sub> analítico ''MC''
! ''θ''<sub>''loc''</sub> numérico ''MC''
|-

!  0<sup>0</sup>
| '''''t'''''<sub>''y''</sub> > 0 y '''''t'''''<sub>''x''</sub> > 0
|  0<sup>''o''</sup>
|  0 . 00<sup>''o''</sup>
|-

!  0<sup>0</sup>
| '''''t'''''<sub>''y''</sub> > 0 y '''''t'''''<sub>''x''</sub> < 0
|  45<sup>''o''</sup>
|  44 . 95<sup>''o''</sup>
|-

!  30<sup>0</sup>
| '''''t'''''<sub>''y''</sub> > 0 y '''''t'''''<sub>''x''</sub> > 0
|  0<sup>''o''</sup>
|  0 . 00<sup>''o''</sup>
|-

!  30<sup>0</sup>
| '''''t'''''<sub>''y''</sub> > 0 y '''''t'''''<sub>''x''</sub> < 0
|  30<sup>''o''</sup>
|  30 . 84<sup>''o''</sup>
|}

La [[#tbl0030|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.

<span id='tbl0030'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 6.

Ángulo de localización ''θ''<sub>''loc''</sub> para el caso uniaxial puro '''''t'''''<sub>''y''</sub> ≠ 0 y '''''t'''''<sub>''x''</sub> = 0 en deformación plana para valores distintos de coeficiente de Poisson

|-

! ''ϕ''
! ''ν''
! ''θ''<sub>''loc''</sub> analítico ''DP''
! ''θ''<sub>''loc''</sub> numérico ''DP''
! ''θ''<sub>''loc''</sub>
! ''θ''<sub>''loc''</sub> numérico ''MC''
|-

!  0<sup>0</sup>
| 0.0
|  45<sup>''o''</sup>
|  44 . 94<sup>''o''</sup>
|  45<sup>''o''</sup>
|  44 . 94<sup>''o''</sup>
|-

!  0<sup>0</sup>
| 0.15
|  45<sup>''o''</sup>
|  44 . 12<sup>''o''</sup>
|  45<sup>''o''</sup>
|  44 . 94<sup>''o''</sup>
|-

!  30<sup>0</sup>
| 0.0
|  28 . 15<sup>''o''</sup>
|  29 . 65<sup>''O''</sup>
|  30<sup>''o''</sup>
|  29 . 65<sup>''o''</sup>
|-

!  30<sup>0</sup>
| 0.15
|  28 . 15<sup>''o''</sup>
|  29 . 65<sup>''o''</sup>
|  30<sup>''o''</sup>
|  29 . 65<sup>''o''</sup>
|}

Finalmente, las figuras [[#fig0055|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.

<span id='fig0055'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr11.jpg|center|px|Curvas Desplazamiento-Reacción.]]


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

Figura 11.

Curvas Desplazamiento-Reacción.

</span>
|}

===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.95''m'', radio externo ''r'' = 0.50''m'' y espesor ''t'' = 0.05''m''; 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 [[#tbl0010|tabla 2]] ( [[#bib0035|[7]]]; [[#bib0050|[10]]] ;  [[#bib0170|[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 [[#fig0060|figura 12]]. En las simulaciones se usan las propiedades materiales: densidad ''ρ'' = 100 ''Kg''/''m''<sup>3</sup>, módulo de Young ''E'' = 10<sup>5</sup>''Pa'', coeficiente de Poisson ''ν'' = 0.30, límite elástico ''f''<sub>''y''</sub> = ''c''<sub>0</sub> = 400 ''Pa  '' y energía de fractura <math display="inline">G_f=5</math> ''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''.

<span id='fig0060'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr12.jpg|center|339px|Geometría y malla empleada.]]


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

Figura 12.

Geometría y malla empleada.

</span>
|}

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 [[#fig0065|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<sup>''o''</sup> con el plano horizontal, con la dirección de la tensión principal menor, ([[#bib0035|[7]]]; [[#bib0050|[10]]] ;  [[#bib0170|[34]]]). La solución numérica coincide exactamente con la solución analítica.

<span id='fig0065'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr13.jpg|center|501px|Deformaciones plásticas equivalentes en el modelo de VM en tracción pura. Ángulo ...]]


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

Figura 13.

Deformaciones plásticas equivalentes en el modelo de VM en tracción pura. Ángulo de localización teórico: <math display="inline">{\theta }_{loc}^{teo}=\pm 45.0^o</math> con el eje vertical. Ángulo de localización numérico: <math display="inline">{\theta }_{loc}^{num}=\pm 45.0^o</math>.

</span>
|}

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<sup>''o''</sup>. 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<sup>''o''</sup> con el plano horizontal. Según los resultados analíticos en la [[#tbl0010|tabla 2]] ([[#bib0035|[7]]]; [[#bib0050|[10]]] ;  [[#bib0170|[34]]]), las líneas de deslizamiento se forman a ±22 . 5<sup>''o''</sup> con la dirección principal menor, es decir, a 22 . 5<sup>''o''</sup> y 67 . 5<sup>''o''</sup> con el plano horizontal. Ambas soluciones alternativas se muestran en las figuras [[#fig0070|14]] y [[#fig0075|15]], respectivamente. Para obtener una u otra, se perturba ligeramente el problema de torsión pura con una tracción/compresión longitudinal, respectivamente.

<span id='fig0070'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr14.jpg|center|503px|Deformaciones plásticas equivalentes en el modelo de MC con ϕ=45o en torsión con ...]]


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

Figura 14.

Deformaciones plásticas equivalentes en el modelo de MC con ''ϕ'' = 45<sup>''o''</sup> en torsión con perturbación de tracción. Ángulo de localización teórico: <math display="inline">{\theta }_{loc}^{teo}=22.5^o</math> con el plano horizontal. Ángulo de localización numérico: <math display="inline">{\theta }_{loc}^{num}=23.{62}^o</math>.

</span>
|}

<span id='fig0075'></span>

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


[[Image:draft_Content_128075237-1-s2.0-S0213131516300244-gr15.jpg|center|503px|Deformaciones plásticas equivalentes en el modelo de MC con ϕ=45o en torsión ...]]


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

Figura 15.

Deformaciones plásticas equivalentes en el modelo de MC con ''ϕ'' = 45<sup>''o''</sup> en torsión pura. Ángulo de localización teórico: <math display="inline">{\theta }_{loc}^{teo}=67.5^o</math> con el plano horizontal. Ángulo de localización numérico: <math display="inline">{\theta }_{loc}^{num}=67.{53}^o</math>.

</span>
|}

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 2''D'' y 3''D'' 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''<sup>''o''</sup> 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==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[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</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] Ivo Babuska; Error-bounds for finite element method; Numerische Mathematik, 16 (4) (1971), pp. 322–333</li>
<li><span id='bib0015'></span>
[[#bib0015|[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</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] M. Cervera, M. Chiumenti; Size effect and localization in {J2} plasticity; International Journal of Solids and Structures, 46 (17) (2009), pp. 3301–3312</li>
<li><span id='bib0025'></span>
[[#bib0025|[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</li>
<li><span id='bib0030'></span>
[[#bib0030|[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</li>
<li><span id='bib0035'></span>
[[#bib0035|[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</li>
<li><span id='bib0040'></span>
[[#bib0040|[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</li>
<li><span id='bib0045'></span>
[[#bib0045|[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</li>
<li><span id='bib0050'></span>
[[#bib0050|[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</li>
<li><span id='bib0055'></span>
[[#bib0055|[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</li>
<li><span id='bib0060'></span>
[[#bib0060|[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</li>
<li><span id='bib0065'></span>
[[#bib0065|[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</li>
<li><span id='bib0070'></span>
[[#bib0070|[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</li>
<li><span id='bib0075'></span>
[[#bib0075|[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</li>
<li><span id='bib0080'></span>
[[#bib0080|[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</li>
<li><span id='bib0085'></span>
[[#bib0085|[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</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] F. Brezzi, M. Fortin, D. Marini; Mixed finite element methods; Springer (1991)</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] GiD. The personal pre and post processor. http://gid.cimne.upc.es. 2009.</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] T.J.R. Hughes, G. Scovazzi, L.P. Franca; Multiscale and Stabilized Methods; John Wiley & Sons, Ltd (2004)</li>
<li><span id='bib0105'></span>
[[#bib0105|[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</li>
<li><span id='bib0110'></span>
[[#bib0110|[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</li>
<li><span id='bib0115'></span>
[[#bib0115|[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</li>
<li><span id='bib0120'></span>
[[#bib0120|[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&rsquo;í (2015)</li>
<li><span id='bib0125'></span>
[[#bib0125|[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</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] E.d.S. Neto, D.R.J. Peric; Computational Methods for Plasticity; Wiley (2008)</li>
<li><span id='bib0135'></span>
[[#bib0135|[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</li>
<li><span id='bib0140'></span>
[[#bib0140|[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</li>
<li><span id='bib0145'></span>
[[#bib0145|[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</li>
<li><span id='bib0150'></span>
[[#bib0150|[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</li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] J.C. Simo, T.J.R. Hughes; Computational Inelasticity; Springer-Verlag, New-York (1998)</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] P. Steinmann, K. Willam; Finite element analysis of elastoplastic discontinuities; Journal of Engineering Mechanics, 120 (11) (1994), pp. 2428–2442</li>
<li><span id='bib0165'></span>
[[#bib0165|[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&rsquo;í, 21 (4) (2005)</li>
<li><span id='bib0170'></span>
[[#bib0170|[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</li>
</ol>