Stabilized mixed explicit finite element formulation for compressible and nearly-incompressible solids

==Resumen==

El presente estudio presenta una formulación mixta de elementos finitos capaz de abordar problemas quasi-incompresibles en forma explícita. Esta formulación se aplica a elementos con interpolaciones independientes e iguales de desplazamientos y deformaciones, estabilizadas mediante subescalas variacionales (VMS). Como continuación del estudio presentado en la referencia [[#bib0115|[23]]] , en la que se introdujo la sub-escala de las deformaciones, en este trabajo se incluyen los efectos de la sub-escala de los desplazamientos, con el fin de estabilizar el campo de las presiones. La formulación evita la condición de Ladyzhenskaya-Babuska-Brezzi (LBB) y sólo requiere la resolución de un sistema diagonal de ecuaciones. En este artículo se tratan también los principales aspectos de implementación. Finalmente, ejemplos de validación numérica muestran el comportamiento de estos elementos en comparación con la formulación irreducible.

==Abstract==

This study presents a mixed finite element formulation able to address nearly-incompressible problems explicitly. This formulation is applied to elements with independent and equal interpolation of displacements and strains, stabilized by variational subscales (VMS). As a continuation of the study presented in reference [[#bib0115|[23]]] , in which the strains sub-scale was introduced, in this work the effects of sub-scale displacements are included, in order to stabilize the pressure field. The formulation avoids the Ladyzhenskaya-Babuska-Brezzi (LBB) condition and only requires the solution of a diagonal system of equations. The main aspects of implementation are also discussed. Finally, numerical examples validate the behaviour of these elements compared with the irreductible formulation.  

==Palabras clave==

Estabilización ; Formulación irreducible ; Formulación mixta explícita

==Keywords==

Stabilization ; Irreductible formulation ; Explicit mixed formulation

==1. Introducción==

En mecánica computacional de sólidos y estructuras es habitual elegir los desplazamientos como incógnitas y calcular las deformaciones por diferenciación del campo de los desplazamientos. Dado que no es posible realizar una reducción adicional en la elección del campo incógnita, los enfoques basados en la formulación estándar de desplazamientos se conocen como ''formulaciones irreducibles'' . Estas formulaciones son muy efectivas en el tratamiento de una amplia gama de problemas prácticos; sin embargo, en ciertos casos, conducen a un comportamiento rígido o dependiente de la malla. Por ejemplo, el rendimiento de los elementos estándares de bajo orden es pobre en condiciones de elasticidad y plasticidad quasi-incompresibles, mostrando un ''bloqueo volumétrico''  que se manifiesta como un comportamiento exageradamente rígido. Existe una vasta literatura proponiendo soluciones para este tipo de problemas, con las primeras propuestas basadas en técnicas de integración reducida o en el uso de deformaciones asumidas, véase, por ejemplo, el método B-bar  [[#bib0090|[18]]]  ;  [[#bib0095|[19]]] .      

A partir los años 90 se emplearon enfoques mixtos basados en desplazamiento-presión ('''''u''''' , ''p'' ) para la solución de problemas incompresibles  [[#bib0030|[6]]] ; [[#bib0085|[17]]] ; [[#bib0120|[24]]] ; [[#bib0130|[26]]] ; [[#bib0140|[28]]]  ;  [[#bib0145|[29]]] . Tales técnicas demostraron también ser aplicables en problemas compresibles (véase, por ejemplo [[#bib0135|[27]]] ). La razón para utilizar formulaciones mixtas es, fundamentalmente, obtener un mejor orden de convergencia en las variables de interés y, especialmente, en las tensiones.      

En este trabajo se emplea una formulación mixta dentro de un esquema explícito de integración temporal y estabilizada mediante el Método de Multi-Escala Variacional (VMS) [[#bib0045|[9]]] ; [[#bib0100|[20]]] ; [[#bib0105|[21]]] ; [[#bib0110|[22]]]  ;  [[#bib0125|[25]]] , a fin de abordar satisfactoriamente problemas de mecánica de sólidos compresibles y quasi-incompresibles, tanto en análisis quasi-estáticos como dinámicos. Aunque no es posible resolver el límite incompresible mediante el empleo de un esquema puramente explícito, en este trabajo se incluyen las subescalas de los desplazamientos a fin de estabilizar el campo de las presiones y proporcionar una solución quasi-incompresible.      

La organización del presente artículo es la siguiente. En la sección [[#sec0010|2]]  se presentan el planteamiento mixto estándar, su correspondiente planteamiento en multi-escalas, los métodos de estabilización y los aspectos relevantes de la implementación numérica de la formulación propuesta. En la sección [[#sec0045|3]]  se comprueba el buen comportamiento de los elementos mixtos explícitos mediante diversas verificaciones numéricas y se compara con los resultados de formulación estándar.      

<span id='sec0010'></span>
==2. Formulación==

===2.1. Formulación mixta===

En la mecánica de sólidos deformables, las deformaciones '''''ɛ'''''  pueden ser consideradas como variables independientes, adicional al campo de los desplazamientos '''''u''''' . Los grados de libertad para un problema tridimensional se incrementan de 3, '''''u'''''  = {''u''<sub>''x''</sub> , ''u''<sub>''y''</sub> , ''u''<sub>''z''</sub> }, a 9, incluyendo los componentes del tensor simétrico de las deformaciones, '''''ɛ'''''  = {''ɛ''<sub>''xx''</sub> , ''ɛ''<sub>''yy''</sub> , ''ɛ''<sub>''zz''</sub> , ''ɛ''<sub>''xy''</sub> , ''ɛ''<sub>''yz''</sub> , ''ɛ''<sub>''xz''</sub> }.      

En tal caso, la formulación ''fuerte''  del problema puede escribirse como: hallar el campo de los desplazamientos '''''u''''' , sus respectivas derivadas temporales (velocidad <math display="inline">\dot{\boldsymbol{u}}</math>  y aceleración <math display="inline">\ddot{\boldsymbol{u}}</math> ) y el campo de deformaciones '''''ɛ''''' , dadas unas fuerzas másicas '''''b''''' , tales que:

<span id='eq0005'></span>
{| 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)
|}

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\varepsilon }-{\nabla }^s\boldsymbol{\mbox{u}}=</math><math>\boldsymbol{0}\quad en\quad \Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

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, expresado como:

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

donde '''''C'''''  representa el tensor constitutivo secante. Además de las ecuaciones [[#eq0005|(1)]]  y [[#eq0010|(2)]] , que se deben satisfacer para cualquier tiempo <math display="inline">t\in I</math>  del intervalo de tiempo de interés <math display="inline">I=(0,T)</math> , el problema está sujeto a unas apropiadas condiciones de contorno, conocidas como condición de Dirichlet y condición de Von Neumann especificadas en el contorno del dominio ∂Ω = ∂Ω<sub>''u''</sub>  ∪ ∂Ω<sub>''σ''</sub>  y ∂Ω<sub>''u''</sub> ∩ ∂Ω<sub>''σ''</sub>  = ∅. Las variables del problema '''''ɛ''''' , '''''u'''''  y sus derivadas, <math display="inline">\dot{\boldsymbol{u}}</math>  están sujetas a unas condiciones iniciales en ''t''  = ''t''<sub>0</sub> , es decir; <math display="inline">\boldsymbol{\varepsilon }\vert _{t=t_0}={\boldsymbol{\varepsilon }}_0=</math><math>{\nabla }^s\boldsymbol{u}_0</math> , <math display="inline">\boldsymbol{u}\vert _{t=t_0}=\boldsymbol{u}_0</math>  y <math display="inline">\dot{\boldsymbol{u}}\vert _{t=t_0}=\boldsymbol{v}_0</math> .      

Después de multiplicar las ecuaciones [[#eq0005|(1)]]  y [[#eq0010|(2)]]  por sus respectivas funciones de peso e integrando por partes la ecuación [[#eq0005|(1)]] , la formulación variacional del problema es:

<span id='eq0020'></span>
{| 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 \quad \forall \boldsymbol{\omega }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

<span id='eq0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }\boldsymbol{\gamma }\cdot \boldsymbol{\varepsilon }\quad d\Omega =</math><math>{\int }_{\Omega }\boldsymbol{\gamma }\cdot {\nabla }^s\boldsymbol{\mbox{u}}\quad d\Omega \quad \quad \forall \boldsymbol{\gamma }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

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 y <math display="inline">\overline{\boldsymbol{t}}=\boldsymbol{\sigma }\boldsymbol{n}</math>  son las tracciones en la dirección normal saliente '''''n'''''  en ∂Ω<sub>''σ''</sub> .      

===2.2. Discretización espacial y temporal===

Para definir la forma discreta del problema, se considera un campo de desplazamientos discretos '''''u'''''  = '''''u'''''<sub>''h''</sub> , sus derivadas temporales discretas (<math display="inline">\dot{\boldsymbol{u}}={\dot{\boldsymbol{u}}}_h</math>  y <math display="inline">\ddot{\boldsymbol{u}}={\ddot{\boldsymbol{u}}}_h</math> ) y un campo de deformaciones discretas '''''ɛ'''''  = '''''ɛ'''''<sub>''h''</sub> . La forma discreta de [[#eq0020|(4)]]  y [[#eq0025|(5)]]  es:

<span id='eq0030'></span>
{| 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 }{\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 \quad \forall {\boldsymbol{\omega }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot {\boldsymbol{\varepsilon }}_h\quad d\Omega =</math><math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot {\nabla }^s\boldsymbol{\mbox{u}}_h\quad d\Omega \quad \quad \forall {\boldsymbol{\gamma }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

donde ahora <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 prueba definidas sobre el espacio de los elementos finitos <math display="inline">V_h</math>  y <math display="inline">T_h</math> , respectivamente. La condición de Ladyzhenskaya-Babuska-Brezzi (''condición inf-sup''  LBB) impone restricciones en las funciones de interpolación para garantizar estabilidad y unicidad de la solución  [[#bib0005|[1]]]  ;  [[#bib0075|[15]]] . De ahí que los espacios de solución para los desplazamientos y las deformaciones no puedan ser escogidos libremente. Lamentablemente, muchos elementos de bajo orden, incluyendo el elemento con igual orden lineal de interpolación para los desplazamientos y las deformaciones, de interés en este trabajo, no cumplen la condición LBB, por lo que no son estables y presentan oscilaciones en el campo de los desplazamientos que contaminan la solución. Afortunadamente, las restricciones de la condición LBB pueden ser evitadas mediante la utilización de apropiadas técnicas de estabilización numérica que pueden proveer la estabilidad necesaria para la convergencia. Los métodos de elementos finitos estabilizados consisten en añadir a la forma discreta términos que son función de los residuos de las ecuaciones [[#eq0030|(6)]]  y [[#eq0035|(7)]]  evaluados dentro de los elementos y que disminuyen con el refinamiento de la malla.      

Además, como las ecuaciones [[#eq0030|(6)]]  y [[#eq0035|(7)]]  son discretas en el espacio pero continuas en el tiempo, se requiere un algoritmo de discretización temporal. En este trabajo se usa el esquema explícito de diferencias centradas.      

===2.3. Método de los elementos finitos estabilizados===

====2.3.1. Multi-escalas variacionales de resolución====

El método de estabilización empleado en este trabajo se basa en el concepto de ''Sub-escalas o Multi-Escalas Variacionales''  (VMS)  [[#bib0105|[21]]] . La idea fundamental es que al definir una malla de elementos finitos quedan establecidas dos ''escalas''  o ''niveles de resolución'' , una que corresponde a la malla y será captado mediante la aproximación de los elementos finitos (escala gruesa) y otra más fina, que corresponde a la parte que la malla no logra captar y que denominaremos simplemente sub-escala. La solución del problema continuo contiene componentes de ambas escalas. De acuerdo con este concepto, se deben aproximar de alguna manera los efectos de la sub-escala con el objetivo de mejorar las propiedades de estabilidad de la formulación. En base a este concepto, el campo de los desplazamientos y deformaciones se aproximan como:

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

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

Las velocidades <math display="inline">\dot{\boldsymbol{u}}</math>  y las aceleraciones <math display="inline">\ddot{\boldsymbol{u}}</math>  se obtienen por diferenciación temporal de los desplazamientos, o sea:

{| 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;" | ( 10)
|}

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

donde <math display="inline">{\boldsymbol{\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 mientras que <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.      

Al considerar las correspondientes funciones de prueba en la sub-escala <math display="inline">\tilde{\boldsymbol{\omega }}\in \tilde{V}</math>  y <math display="inline">\tilde{\boldsymbol{\gamma }}\in \tilde{T}</math> , nótese que la aproximación de la solución en el campo de deformaciones y desplazamientos se extiende ahora a <math display="inline">T\simeq T_h\oplus \tilde{T}</math>  y <math display="inline">V\simeq V_h\oplus \tilde{V}</math>  respectivamente. Cada método particular de elementos finitos estabilizados se define de acuerdo a cómo se escogen los espacios <math display="inline">\tilde{V}</math>  y <math display="inline">\tilde{T}</math> . El método clásico de Galerkin de obtiene cuando <math display="inline">\tilde{V}=\tilde{T}=\lbrace 0\rbrace </math> .      

Por otra parte, el tensor de tensiones en la ecuación [[#eq0015|(3)]]  toma la forma:

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

en el cual '''''σ'''''<sub>''h''</sub>  = '''''C''''' ('''''ɛ'''''<sub>''h''</sub> ) : '''''ɛ'''''<sub>''h''</sub>  y <math display="inline">\tilde{\boldsymbol{\sigma }}=\boldsymbol{C}({\boldsymbol{\varepsilon }}_h):\tilde{\boldsymbol{\varepsilon }}</math>  son respectivamente las tensiones ''continuas''  y las tensiones de la escala fina. Obsérvese que se ha empleado la aproximación '''''C''''' ('''''ɛ''''' ) ≈ '''''C''''' ('''''ɛ'''''<sub>''h''</sub> ). Esta aproximación y la correspondiente descomposición aditiva del tensor de tensiones en la ecuación [[#eq0060|(12)]]  están justificados en  [[#bib0015|[3]]]  ;  [[#bib0020|[4]]] .      

Introduciendo las ecuaciones [[#eq0040|(8)]] , [[#eq0045|(9)]] , [[#eq0055|(11)]]  y [[#eq0060|(12)]]  en [[#eq0020|(4)]]  y [[#eq0025|(5)]] , el problema se lee ahora como:

<span id='eq0065'></span>
{| 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>\boldsymbol{F}_{ext}({\boldsymbol{\omega }}_h)\quad \forall {\boldsymbol{\omega }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

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

<span id='eq0075'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot {\boldsymbol{\varepsilon }}_h\quad d\Omega +</math><math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot \tilde{\boldsymbol{\varepsilon }}\quad d\Omega =</math><math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot {\nabla }^s\boldsymbol{\mbox{u}}_h\quad d\Omega -</math><math>{\int }_{\Omega }\left(\nabla \cdot {\boldsymbol{\gamma }}_h\right)\cdot \tilde{\boldsymbol{\mbox{u}}}\quad d\Omega \quad \forall {\boldsymbol{\gamma }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

<span id='eq0080'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot {\boldsymbol{\varepsilon }}_h\quad d\Omega +</math><math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot \tilde{\boldsymbol{\varepsilon }}\quad d\Omega =</math><math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot {\nabla }^s\boldsymbol{\mbox{u}}_h\quad d\Omega +</math><math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot {\nabla }^s\tilde{\boldsymbol{\mbox{u}}}\quad d\Omega \quad \forall \tilde{\boldsymbol{\gamma }}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

donde en la ecuación [[#eq0075|(15)]]  se ha empleado la integración por partes y se ha despreciado la integral sobre contorno [[#bib0015|[3]]]  ;  [[#bib0045|[9]]] . Las fuerzas externas '''''F'''''<sub>''ext''</sub> ('''''ω'''''<sub>''h''</sub> ) se escriben de la forma:

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

Debido a la aproximación empleada en [[#eq0040|(8)]]  y [[#eq0045|(9)]] , la ecuación [[#eq0005|(1)]]  se desdobla en las dos ecuaciones [[#eq0065|(13)]]  y [[#eq0070|(14)]] , cada una de ellas relacionadas con la escala correspondiente. Asimismo, la ecuación [[#eq0010|(2)]]  se desdobla en las ecuaciones [[#eq0075|(15)]]  y [[#eq0080|(16)]] .      

Dado que las ecuaciones [[#eq0070|(14)]]  y [[#eq0080|(16)]]  no se resuelven dentro del espacio de los elementos finitos, se requiere una aproximación de las mismas. Obsérvese que no son necesarios los valores punto a punto, sino sólo el valor de la integral correspondiente al término de estabilización en [[#eq0065|(13)]]  y [[#eq0075|(15)]][[#bib0015|[3]]] . El objetivo de estas aproximaciones es obtener una estimación variacional de los efectos de las sub-escalas <math display="inline">\tilde{\boldsymbol{u}}</math>  y <math display="inline">\tilde{\boldsymbol{\varepsilon }}</math>  en las ecuaciones de balance.      

Los residuos correspondientes a las ecuaciones en la sub-escala de los elementos finitos son:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{r}_{\varepsilon ,h}={\nabla }^s\boldsymbol{u}_h-</math><math>{\boldsymbol{\varepsilon }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

<span id='eq0095'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{r}_{u,h}=\rho {\ddot{\boldsymbol{u}}}_h-</math><math>\nabla \cdot {\boldsymbol{\sigma }}_h-\boldsymbol{b}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

Esto permite escribir las ecuaciones [[#eq0070|(14)]]  y [[#eq0080|(16)]]  como:

<span id='eq0100'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }\tilde{\boldsymbol{\omega }}\cdot \left(\nabla \cdot \left(\boldsymbol{C}:\tilde{\boldsymbol{\varepsilon }}\right)\right)\quad d\Omega -</math><math>{\int }_{\Omega }\rho \tilde{\boldsymbol{\omega }}\cdot \ddot{\tilde{\boldsymbol{u}}}\quad d\Omega =</math><math>{\int }_{\Omega }\tilde{\boldsymbol{\omega }}\cdot \boldsymbol{r}_{u,h}\quad d\Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

<span id='eq0105'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot \tilde{\boldsymbol{\varepsilon }}\quad d\Omega -</math><math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot {\nabla }^s\tilde{\boldsymbol{u}}\quad d\Omega =</math><math>{\int }_{\Omega }\tilde{\boldsymbol{\gamma }}\cdot \boldsymbol{r}_{\varepsilon ,h}\quad d\Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

Los términos de la derecha en las ecuaciones [[#eq0100|(20)]]  y [[#eq0105|(21)]]  son las proyecciones <math display="inline">P(\bullet )</math>  de los residuos sobre el espacio de las sub-escalas. Por tanto, las ecuaciones [[#eq0100|(20)]]  y [[#eq0105|(21)]]  relacionan las sub-escalas <math display="inline">\tilde{\boldsymbol{\varepsilon }}</math>  y <math display="inline">\tilde{\boldsymbol{u}}</math>  con los residuos de la escala de elementos finitos, representados por '''''r'''''<sub>''ɛ'' ,''h''</sub>  y '''''r'''''<sub>''u'' ,''h''</sub> , respectivamente.      

Existen diversas posibilidades para definir la aproximación a las sub-escalas. En principio el espacio de la escala fina podría ser cualquier espacio complementario al espacio de elementos finitos. Tómese en cuenta que en la ecuación [[#eq0100|(20)]]  aparece la segunda derivada temporal de la sub-escala <math display="inline">\tilde{\boldsymbol{u}}</math> , por lo que ésta necesita discretizarse en el tiempo [[#bib0060|[12]]] . Es necesario, en primer lugar, seguir la evolución en el tiempo de la sub-escala y, en segundo lugar, el método de integración temporal de la sub-escalas tiene que ser compatible con el esquema explícito empleado en el espacio de los elementos finitos, por lo que también se requiere resolver en forma explícita la sub-escala de los desplazamientos <math display="inline">\tilde{\boldsymbol{u}}</math> . La subescala de las deformaciones <math display="inline">\tilde{\boldsymbol{\varepsilon }}</math>  puede ser tratada implícitamente (véase [[#bib0045|[9]]] ).      

Considérese una partición uniforme en el tiempo [0, ''T'' ] de tamaño Δ''t'' , tal que el tiempo ''t''<sup>''n''</sup>  = ''n'' Δ''t'' . Por simplicidad, se supone Δ''t''  constante. Empleando el método de diferencias centrales en la ecuación  [[#eq0100|(20)]] , las aceleraciones en la sub-escala se escriben:

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

Sustituyendo [[#eq0110|(22)]]  en [[#eq0100|(20)]]  y siendo <math display="inline">P_{\boldsymbol{\varepsilon }}</math>  y <math display="inline">P_\boldsymbol{u}</math>  las proyecciones sobre <math display="inline">\tilde{T}</math>  y <math display="inline">\tilde{V}</math> , respectivamente, el problema [[#eq0100|(20)]]  y [[#eq0105|(21)]]  se lee ahora como: dados <math display="inline">{\tilde{\boldsymbol{u}}}^n</math>  y <math display="inline">{\tilde{\boldsymbol{u}}}^{n-1}</math> , hallar <math display="inline">{\tilde{\boldsymbol{u}}}^{n+1}</math>  y <math display="inline">{\tilde{\boldsymbol{\varepsilon }}}^{n+1}</math>  tal que [[#bib0045|[9]]]  ;  [[#bib0060|[12]]] :

<span id='eq0115'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>P_u\left(\nabla \cdot \boldsymbol{C}({\boldsymbol{\varepsilon }}_h):{\tilde{\boldsymbol{\varepsilon }}}^{n+1}-\right. </math><math>\left. \frac{\rho }{\Delta t^2}{\tilde{\boldsymbol{u}}}^{n+1}\right)=</math><math>P_u\left(\frac{\rho }{\Delta t^2}(2{\tilde{\boldsymbol{u}}}^n-\right. </math><math>\left. {\tilde{\boldsymbol{u}}}^{n-1})+\boldsymbol{r}_{u,h}^n\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)
|}

<span id='eq0120'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>P_{\boldsymbol{\varepsilon }}\left({\tilde{\boldsymbol{\varepsilon }}}^{n+1}-\right. </math><math>\left. {\nabla }^s{\tilde{\boldsymbol{u}}}^{n+1}\right)=</math><math>P_{\boldsymbol{\varepsilon }}\left(\boldsymbol{r}_{\varepsilon ,h}^{n+1}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)
|}

Se requiere resolver de forma aproximada las ecuaciones [[#eq0115|(23)]]  y [[#eq0120|(24)]]  para obtener los valores de <math display="inline">{\tilde{\boldsymbol{\varepsilon }}}^{n+1}</math><math display="inline">{\tilde{\boldsymbol{u}}}^{n+1}</math> . Aplicando el mismo procedimiento que en [[#bib0055|[11]]] , se pueden aproximar las sub-escalas mediante un análisis de Fourier del problema. Así, <math display="inline">{\tilde{\boldsymbol{\varepsilon }}}^{n+1}</math>  y <math display="inline">{\tilde{\boldsymbol{u}}}^{n+1}</math>  pueden aproximarse dentro de cada elemento Ω<sub>''e''</sub>  como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\tilde{\boldsymbol{u}}}}^{n+1}\approx {\tau }_{u_t}\frac{\rho }{\Delta t^2}\left(2{\tilde{\boldsymbol{u}}}^n-\right. </math><math>\left. {\tilde{\boldsymbol{u}}}^{n-1}\right)+{\tau }_{u_t}P_u\left(-\right. </math><math>\left. \boldsymbol{r}_{u,h}^n\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tilde{\boldsymbol{\varepsilon }}}^{n+1}\approx {\tau }_{\varepsilon }P_{\varepsilon }\left(\boldsymbol{r}_{\varepsilon ,h}^{n+1}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)
|}

donde <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}+\frac{1}{{\tau }_{u_s}}\right)}^{-1}\quad y\quad {\tau }_{\varepsilon }=</math><math>c_{\varepsilon }\frac{h}{L_0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)
|}

donde <math display="inline">{\tau }_{u_s}=c_u{hL}_0/\mu </math>  es el parámetro de estabilización estático, ''c''<sub>''u''</sub>  > 0 y ''c''<sub>''ɛ''</sub>  > 0 son constantes algorítmicas adimensionales, ''L''<sub>0</sub>  es la longitud característica del problema el cual se toma como el diámetro del dominio computacional Ω y ''μ''  > 0 es el módulo de rigidez secante al corte. La experiencia numérica en análisis estáticos muestra que los valores de ''c''<sub>''ɛ''</sub>  y ''c''<sub>''u''</sub>  se pueden escoger en el rango de [0.01, 1.0].      

Para definir completamente el método y calcular las sub-escalas, se requiere seleccionar adecuadamente las proyecciones <math display="inline">P_\boldsymbol{u}</math>  y <math display="inline">P_{\boldsymbol{\varepsilon }}</math> . La elección más simple para las proyecciones es tomar las mismas como la identidad <math display="inline">I</math> . Este procedimiento se conoce como ''Método de las Sub-Escalas Algebraicas''  (Algebraic Sub-Scales (ASGS))  [[#bib0040|[8]]] . Codina propuso para el espacio de las sub-escalas <math display="inline">\tilde{V}</math>  y <math display="inline">\tilde{T}</math>  el espacio ortogonal al espacio de los elementos finitos <math display="inline">V_h</math>  y <math display="inline">T_h</math> . Esta definición da origen a una formulación denominada ''Método de las Sub-escalas Ortogonales''  (Orthogonal Sub-Scales (OSS))  [[#bib0040|[8]]] ; [[#bib0045|[9]]]  ;  [[#bib0050|[10]]] . De acuerdo con este planteamiento, el espacio de las deformaciones y desplazamientos se aproxima como <math display="inline">T\simeq T_h\oplus T_h^{\perp }</math>  y <math display="inline">V\simeq V_h\oplus V_h^{\perp }</math>  respectivamente. En este trabajo se ha optado por este método. Por practicidad se introduce el operador de proyección ortogonal <math display="inline">P_h^{\perp }(\bullet )=(\bullet )-P_h(\bullet )</math> . Además se supone que el campo de las fuerzas másicas '''''b'''''  puede ser descrito completamente en el espacio de <math display="inline">V_h</math> , se tiene que <math display="inline">P_h^{\perp }(\boldsymbol{b})=0</math>[[#bib0010|[2]]] ; [[#bib0025|[5]]]  ;  [[#bib0150|[30]]] . Por consiguiente, se pueden tomar <math display="inline">\tilde{\boldsymbol{u}}</math>  y <math display="inline">\tilde{\boldsymbol{\varepsilon }}</math>  como:

<span id='eq0140'></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{\tilde{\boldsymbol{u}}}^n-\right. </math><math>\left. {\tilde{\boldsymbol{u}}}^{n-1}\right)+{\tau }_{u_t}\left(\nabla \cdot {\boldsymbol{\sigma }}_h^n-\right. </math><math>\left. P_h\left(\nabla \cdot {\boldsymbol{\sigma }}_h^n\right)\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)
|}

<span id='eq0145'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tilde{\boldsymbol{\varepsilon }}}^{n+1}\simeq {\tau }_{\varepsilon }\left({\nabla }^s\boldsymbol{\mbox{u}}_h^{n+1}-\right. </math><math>\left. P_h\left({\nabla }^s\boldsymbol{\mbox{u}}_h^{n+1}\right)\right)=</math><math>{\tau }_{\varepsilon }\left({\nabla }^s\boldsymbol{\mbox{u}}_h^{n+1}-\right. </math><math>\left. {\breve{\boldsymbol{\varepsilon }}}_h^{n+1}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)
|}

siendo <math display="inline">{\breve{\boldsymbol{\varepsilon }}}_h^{n+1}=P_h\left({\nabla }^s\boldsymbol{\mbox{u}}_h^{n+1}\right)</math> .      

El residuo en la ecuación [[#eq0095|(19)]]  puede escribirse en función de la parte volumétrica y desviadora del tensor de tensiones, es decir, ∇ · '''''σ'''''<sub>''h''</sub>  = ∇ · '''''S'''''<sub>''h''</sub>  + ∇ ''p''<sub>''h''</sub> , siendo '''''S'''''  el tensor desviador. En problemas quasi-incompresibles o incompresibles, sólo el gradiente de la presión ∇''p''<sub>''h''</sub>  necesita incluirse en el residuo para asegurar la estabilidad en el campo de las presiones. Esto permite estabilizar la presión minimizando la magnitud del término de estabilización. El hecho de no considerar una parte del residuo introduce un cierto error de consistencia, pero en el método de OSS este error es de orden óptimo y no altera la velocidad de convergencia del esquema.      

Por otra parte, se puede introducir un efecto simple de amortiguamiento numérico a fin de considerar un esquema disipativo y eliminar las frecuencias espurias en la sub-escala de los desplazamientos, que pueden afectar la estabilidad del sistema. En tal sentido, despreciando el término desviador en la ecuación [[#eq0140|(28)]]  e introduciendo un coeficiente de disipación <math display="inline">\tilde{\xi }\in [0,1]</math> , se obtiene una versión modificada de OSS [[#bib0010|[2]]] . En consecuencia, <math display="inline">\tilde{\boldsymbol{u}}</math>  estará ahora dado por:

<span id='eq0150'></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;" | ( 30)
|}

===2.4. Ecuación discreta y estabilizada de la ecuación de movimiento. Formulación mixta estabilizada explícita===

Para definir la ecuación discreta y estabilizada de la ecuación de movimiento [[#eq0065|(13)]] , ésta tiene que ser discretizada en el tiempo. En la referencia [[#bib0115|[23]]]  se discuten las ventajas de la formulación mixta explícita en conjunto con un estudio numérico de estabilidad y convergencia. Naturalmente, este método mixto-explícito es condicionalmente estable, pero su precisión es superior a su contraparte irreducible. Al igual que la formulación irreducible, los desplazamientos '''''u'''''  y deformaciones '''''ɛ'''''  en el tiempo ''t''<sup>''n'' +1                            </sup>  se computan a partir de la ecuación de equilibrio dinámico evaluada en el tiempo ''t''<sup>''n''</sup> . En consecuencia, la ecuación [[#eq0065|(13)]]  se expresa en el tiempo ''t''<sup>''n''</sup>  como:

<span id='eq0155'></span>
{| 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^n\quad d\Omega +</math><math>{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:{\tilde{\boldsymbol{\sigma }}}^n\quad d\Omega +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho \left({\ddot{\boldsymbol{u}}}_h^n+\right. </math><math>\left. {\ddot{\tilde{\boldsymbol{u}}}}^n\right)\quad d\Omega =</math><math>\boldsymbol{F}_{ext}^n({\boldsymbol{\omega }}_h)\quad \forall {\boldsymbol{\omega }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)
|}

Al emplear el método de estabilización OSS, el término <math display="inline">{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho {\ddot{\tilde{\boldsymbol{u}}}}^n\quad d\Omega =</math><math>\boldsymbol{0}</math>  por la condición de ortogonalidad. Sustituyendo la ecuación [[#eq0145|(29)]]  en [[#eq0155|(31)]] , la ecuación discreta y estabilizada de equilibrio dinámico en ''t''  = ''t''<sup>''n''</sup>  se escribe finalmente como:

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

Por otra parte, sustituyendo la ecuación [[#eq0145|(29)]]  en la ecuación [[#eq0075|(15)]]  evaluada en ''t''<sup>''n'' +1                            </sup>  se obtiene:

<span id='eq0165'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot {\boldsymbol{\varepsilon }}_h^{n+1}\quad d\Omega +</math><math>{\tau }_{\boldsymbol{\varepsilon }}{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot \left({\nabla }^s\boldsymbol{\mbox{u}}_h^{n+1}-\right. </math><math>\left. {\breve{\boldsymbol{\varepsilon }}}_h^{n+1}\right)\quad d\Omega =</math><math>{\int }_{\Omega }{\boldsymbol{\gamma }}_h\cdot {\nabla }^s\boldsymbol{\mbox{u}}_h^{n+1}\quad d\Omega -</math><math>{\int }_{\Omega }\left(\nabla \cdot {\boldsymbol{\gamma }}_h\right)\cdot {\tilde{\boldsymbol{u}}}^{n+1}\quad d\Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)
|}

Por la propiedad de idempotencia del operador de proyección, es decir, <math display="inline">P_h(P_h(\bullet ))=P_h({\left(\bullet \right)}_h)=</math><math>{\left(\bullet \right)}_h</math> , y teniendo en cuenta que las sub-escalas de los desplazamientos en ''t''<sup>''n'' +1                            </sup>  están dados por la ecuación [[#eq0140|(28)]]  o [[#eq0150|(30)]] , las deformaciones continuas <math display="inline">{\boldsymbol{\varepsilon }}_h^{n+1}</math>  se expresan finalmente:

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

Si no se considera la sub-escala de los desplazamientos, el segundo término del lado derecho de la ecuación [[#eq0170|(34)]]  es nulo, obteniéndose la misma expresión que en [[#bib0115|[23]]] . Por otro lado, dado que el segundo término del lado derecho de la ecuación [[#eq0170|(34)]]  es pequeño frente a <math display="inline">{\breve{\boldsymbol{\varepsilon }}}_h^{n+1}</math> , se puede tomar <math display="inline">{\breve{\boldsymbol{\varepsilon }}}_h\approx {\boldsymbol{\varepsilon }}_h</math> . En tal caso, la ecuación de movimiento estabilizada de equilibrio dinámico en [[#eq0165|(33)]]  se reescribe como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>(1-{\tau }_{\varepsilon }){\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:\boldsymbol{C}({\boldsymbol{\varepsilon }}_h^n):{\boldsymbol{\varepsilon }}_h^n\quad d\Omega +</math><math>{\tau }_{\varepsilon }{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:\boldsymbol{C}({\boldsymbol{\varepsilon }}_h^n):\left({\nabla }^s\boldsymbol{\mbox{u}}_h^n\right)\quad d\Omega +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho {\ddot{\boldsymbol{u}}}_h^n\quad d\Omega =</math><math>\boldsymbol{F}_{ext}^n({\boldsymbol{\omega }}_h)\quad \quad \quad \forall {\boldsymbol{\omega }}_h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)
|}

Reagrupando los términos, se obtiene la expresión:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\varepsilon }}_h^{stab,n}=(1-{\tau }_{\varepsilon }){\boldsymbol{\varepsilon }}_h^n-</math><math>{\tau }_{\varepsilon }{\nabla }^s\boldsymbol{u}_h^n</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)
|}

<span id='eq0185'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:\boldsymbol{C}({\boldsymbol{\varepsilon }}_h^n):{\boldsymbol{\varepsilon }}_h^{stab,n}\quad d\Omega +</math><math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho {\ddot{\boldsymbol{u}}}_h^n\quad d\Omega =</math><math>\boldsymbol{F}_{ext}^n({\boldsymbol{\omega }}_h)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)
|}

donde '''''ɛ'''''<sup>''stab''</sup>  son las ''deformaciones estabilizadas'' . El primer término del lado izquierdo de la la ecuación  [[#eq0185|(37)]]  es el vector de fuerzas internas estabilizadas:

<span id='eq0190'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{F}_{int}^{stab,n}({\boldsymbol{\omega }}_h)=</math><math>{\int }_{\Omega }{\nabla }^s{\boldsymbol{\omega }}_h:{\boldsymbol{\sigma }}^{stab,n}\quad d\Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)
|}

siendo '''''σ'''''<sup>''stab''</sup>  = '''''C'''''  : '''''ɛ'''''<sup>''stab''</sup>  las ''tensiones estabilizadas'' . Esto permite reescribir la ecuación  [[#eq0185|(37)]]  como:

<span id='eq0195'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }{\boldsymbol{\omega }}_h\cdot \rho {\ddot{\boldsymbol{u}}}_h^n\quad d\Omega =</math><math>\boldsymbol{F}_{ext}^n({\boldsymbol{\omega }}_h)-\boldsymbol{F}_{int}^{stab,n}({\boldsymbol{\omega }}_h)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)
|}

Si se considera las fuerzas viscosas, la ecuación [[#eq0195|(39)]]  toma el formato matricial de equilibrio dinámico:

<span id='eq0200'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{M}{\ddot{\boldsymbol{u}}}_h^n+\boldsymbol{D}{\dot{\boldsymbol{u}}}_h^n=</math><math>\boldsymbol{F}_{ext}^n-\boldsymbol{F}_{int}^{stab,n}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)
|}

donde '''''M'''''  es la matriz de masa y '''''D'''''  la matriz de amortiguamiento. El proceso de integración temporal de la ecuación de movimiento en un esquema mixto explícito es idéntico al procedimiento realizado en un esquema irreducible, salvo que se requieren fundamentalmente dos pasos adicionales (3 y 4), los cuales se describen a continuación:
* 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 [[#eq0140|(28)]]  o [[#eq0150|(30)]] .                          
* 4. Actualización de las deformaciones <math display="inline">{\boldsymbol{\varepsilon }}_h^{n+1}</math>  con [[#eq0170|(34)]] .                          
* 5. Ir a siguiente paso.

La formulación propuesta puede ser incorporada con facilidad en cualquier código de elementos finitos explícitos agregando simplemente los pasos 3 y 4, que pueden ejecutarse en única rutina de cálculo. Además, en los pasos 3 y 4 pueden aplicarse técnicas de paralización computacional ya que las operaciones involucradas se realizan a nivel de elementos y nodos de la malla.

Nótese que los desplazamientos <math display="inline">\boldsymbol{u}_h^{n+1}</math>  se obtienen a partir de los valores del paso ''t''<sup>''n''</sup>  y más específicamente de <math display="inline">{\boldsymbol{\varepsilon }}_h^{stab,n}</math> , por lo que la ecuación [[#eq0170|(34)]]  es compatible con un esquema explícito. Una vez calculadas los valores de los desplazamientos <math display="inline">\boldsymbol{u}_h^{n+1}</math> , se calculan los valores de los desplazamientos de la sub-escala <math display="inline">{\tilde{\boldsymbol{u}}}^{n+1}</math>  y por último los valores de las deformaciones nodales <math display="inline">{\boldsymbol{\varepsilon }}_h^{n+1}</math>  con la ecuación [[#eq0170|(34)]] . Este esquema sólo requiere la resolución de un sistema diagonal de ecuaciones, frente a las formulaciones mixtas implícitas, que requieren la resolución de mayores sistemas de ecuaciones. A partir de ahora, a esta formulación propuesta la denominaremos ''Método de Elementos Finitos Mixtos Explícitos''  (Mixed Explicit Finite Element Method) o simplemente MEX-FEM, por sus siglas en inglés.      

===2.5. Aspectos de implementación computacional===

La integración temporal explícita es efectiva si las matrices de masa '''''M'''''  y de amortiguamiento '''''D'''''  en [[#eq0200|(40)]]  son diagonales [[#bib0035|[7]]] . Considérese un elemento tetraedro de cuatros nodos empleando las mismas funciones lineales de interpolación para los desplazamientos y deformaciones. En la formulación mixta, las deformaciones no son constantes dentro del elemento, sino que varían linealmente. En este trabajo se usa una cuadratura cerrada Gauss-Lobatto, en la cual los puntos de integración se sitúan sobre los nodos de los elementos finitos, a fin de evaluar correctamente las integrales sobre el dominio de los elementos finitos Ω<sub>''e''</sub> .      

Nótese que que al emplear el método de Galerkin, las funciones de prueba '''''ω'''''<sub>''h''</sub>  para el espacio de los desplazamientos y '''''γ'''''<sub>''h''</sub>  para el espacio de las deformaciones son las funciones de interpolación, convencionalmente escritas en notación matricial como <math display="inline">\boldsymbol{N}_{\alpha }=[\boldsymbol{N}_{\alpha }^1\boldsymbol{N}_{\alpha }^2\boldsymbol{N}_{\alpha }^3\boldsymbol{N}_{\alpha }^4]</math>  ∀''α''  ∈ {'''''u''''' , '''''ɛ''''' }. Los desplazamientos y las deformaciones nodales dentro del elemento se interpolan como:

{| 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;" | ( 41)
|}

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 otra parte, en la ecuación [[#eq0190|(38)]] , <math display="inline">{\nabla }^s{\boldsymbol{\omega }}_h=\boldsymbol{B}_u^T</math>  es 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;" | ( 42)
|}

De igual manera, el operador de divergencia ∇ · '''''γ'''''<sub>''h''</sub>  en [[#eq0170|(34)]]  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> .      

<span id='sec0045'></span>
==3. Verificación. Simulaciones numéricas==

La eficacia de la formulación MEX-FEM se muestra a continuación en una serie de tests numéricos realizados sobre la Membrana de Cook, mostrado en la [[#fig0005|figura 1]] . Se consideran casos de elasticidad compresible y quasi-incompresible, con las constantes elásticas: Módulo de Young ''E''  = 200 Mpa y coeficiente de Poisson ''ν''  = 0.30 y ''ν''  = 0.499. La densidad del material se toma 10 kg/m<sup>3</sup> . El problema se analiza realizando análisis quasi-estático y dinámico bajo la hipótesis de deformación plana (con espesor unitario) y haciendo un análisis completamente tridimensional (espesor 10 mm). Para el caso bidimensional, el tamaño del elemento finito se ha tomado como <math display="inline">h={\left(4/\pi \cdot A\right)}^{\frac{1}{2}}</math>  siendo ''A  ''  el área del elemento finito correspondiente. En el caso tridimensional se ha tomado <math display="inline">h={\left(6/\pi \cdot V\right)}^{\frac{1}{3}}</math>  donde ''V''  el volumen del elemento finito. Como método de estabilización se emplea el OSS con ''c''<sub>''u''</sub>  = 1.0, ''c''<sub>''ɛ''</sub>  = 1.0 y ''L''<sub>0</sub>  = 50 mm. La sub-escala de los desplazamientos <math display="inline">\tilde{\boldsymbol{u}}</math>  se computa empleando la ecuación [[#eq0150|(30)]] . En la integración temporal de la ecuación de movimiento se ha empleado un paso de tiempo Δ''t  ''  constante. Los resultados obtenidos con la formulación propuesta se comparan con la formulación estándar en desplazamientos y la formulación mixta explícita estándar (Mex-Fem estándar <math display="inline">\tilde{\boldsymbol{u}}=\boldsymbol{0}</math> ) [[#bib0115|[23]]] , es decir sin considerar la sub-escala de los desplazamientos. Los análisis se han realizado con el programa de elementos finitos KRATOS  [[#bib0065|[13]]]  ;  [[#bib0070|[14]]] , desarrollado en el Centro Internacional de Métodos Numéricos (CIMNE). Como pre y post-procesador se ha utilizado GiD [[#bib0080|[16]]] , también desarrollado en CIMNE.

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr1.jpg|center|269px|Membrana de Cook. Geometría.]]


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

Figura 1.

''Membrana de Cook'' . Geometría.                  

</span>
|}

En la discretización espacial se emplean 4 mallas de elementos finitos triangulares y tetraédricos mostradas en las figuras [[#fig0010|2]]  y [[#fig0015|3]] , respectivamente. Para eliminar las frecuencias espurias que pueden aparecer en la integración temporal sin degradar la respuesta dinámica se emplea un amortiguamiento de Rayleigh con un coeficiente disipación en la sub-escala <math display="inline">\tilde{\xi }=0.1</math> . Asimismo, con el fin de obtener la respuesta estacionaria del problema quasi-estático se emplea un amortiguamiento de Rayleigh con un coeficiente de amortiguamiento de ''ξ''  = 0.1. La carga impuesta ''F''  de valor unitario se aplica de una manera instantánea y constante en el tiempo.

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr2.jpg|center|498px|Membrana de Cook. Mallas de elementos finitos empleada en el análisis ...]]


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

Figura 2.

Membrana de Cook. Mallas de elementos finitos empleada en el análisis bidimensional. Nn: número de elementos y Nn: número de nodos.

</span>
|}

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr3.jpg|center|492px|Membrana de Cook. Mallas de elementos finitos empleada en el análisis ...]]

Figura 3.

Membrana de Cook. Mallas de elementos finitos empleada en el análisis tridimensional. Nn: número de elementos y Nn: número de nodos.

La evolución del desplazamiento vertical en la esquina superior derecha (punto A) para un análisis 2''D''  de deformación plana se muestra en las figuras  [[#fig0020|4]]  y [[#fig0025|5]]  para valores de Poisson de ''ν''  = 0.3 y ''ν''  = 0.499, casos compresible y quasi-incompresible, respectivamente. Se puede observar que para ''ν''  = 0.3, ambas formulaciones presentan resultados similares en el campo de los desplazamientos y presiones, aunque MEX-FEM es más preciso. Sin embargo, nótese como en el caso quasi-incompresible ''ν''  = 0.499, la formulación propuesta claramente supera a la irreducible. La [[#fig0025|figura 5]] a muestra el pobre comportamiento del elemento triangular estándar en el campo de los desplazamientos debido al efecto de bloqueo volumétrico. El desplazamiento predicho por este elemento, incluso utilizando mallas finas, está muy por debajo del valor correcto. Nótese además que el error de fase en la formulación irreducible es notoriamente mayor que el obtenido con la presente formulación para una misma malla de elementos finitos. La evolución del campo de las presiones para ''ν''  = 0.499 (punto ''B''  de la  [[#fig0005|figura 1]] ) se aprecia en las figuras [[#fig0025|5]] c y [[#fig0025|5]] d, tanto para la formulación irreducible como para la presente formulación. Dado que el valor de la presión es constante dentro del elemento lineal irreducible, el valor nodal se ha obtenido mediante la proyección <math display="inline">L_2</math> , o sea, <math display="inline">p_h=P_h(p)</math> . Obsérvese cómo el efecto del bloqueo volumétrico en la formulación irreducible provoca oscilaciones el campo de las presiones que desvirtúan completamente la solución. En la [[#fig0025|figura 5]] d se puede comprobar que la presente formulación produce resultados satisfactorios y razonablemente precisos en el campo de las presiones.

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr4.jpg|center|px|Membrana de Cook 2D. Formulación Irreducible vs MEX-FEM con ν=0.3.]]


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

Figura 4.

Membrana de Cook 2''D'' . Formulación Irreducible vs MEX-FEM con ''ν''  = 0.3.                  

</span>
|}

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr5.jpg|center|px|Membrana de Cook 2D. Formulación Irreducible vs MEX-FEM con ν=0.499.]]


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

Figura 5.

Membrana de Cook 2''D'' . Formulación Irreducible vs MEX-FEM con ''ν''  = 0.499.                  

</span>
|}

Las curvas de convergencia para el caso estacionario para los valores de coeficiente de Poisson mencionados se muestran en la [[#fig0030|figura 6]] . Nuevamente para valor ''ν''  = 0.3, la formulación irreducible, MEX-FEM con <math display="inline">\tilde{\boldsymbol{u}}=0</math>  y MEX-FEM convergen al valor correcto del desplazamiento (1.843 mm) y presión (1.632 MPa) a media que se refina la malla. Además, puede verse que las formulaciones mixta estándar y la propuesta convergen a la misma velocidad pero más rápido que la formulación irreducible. Por otra parte, cuando ''ν''  = 0.499, se observa la capacidad de MEX-FEM para predecir el valor correcto del desplazamiento vertical (1.554 mm) y la presión (1.872 MPa). Obsérvese también que la formulación mixta explícita estándar, a pesar de dar un resultado satisfactorio del campo de los desplazamientos, presenta inestabilidades en el campo de las presiones que desvirtúan completamente la solución. En el mismo sentido, nótese que el valor del desplazamiento y de la presión predicho por la formulación irreducible están subestimados.

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr6.jpg|center|px|Membrana de Cook 2D. Curvas de convergencia.]]


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

Figura 6.

Membrana de Cook 2''D'' . Curvas de convergencia.                  

</span>
|}

La [[#fig0035|figura 7]]  muestra el campo de las presiones de las formulaciones mencionadas en la Malla ''D''  para un valor de ''ν''  = 0.499. Obsérvese como MEX-FEM es capaz de dar una solución estabilizada del campo de las presiones, no siendo así para el caso irreducible y el caso mixto explícito sin sub-escala de desplazamientos.

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr7.jpg|center|501px|Membrana de Cook 2D. Campo de presiones. Malla D.]]


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

Figura 7.

Membrana de Cook 2''D'' . Campo de presiones. Malla D.                  

</span>
|}

Los resultados para el análisis dinámico tridimensional quasi-incompresible de la membrana de Cook se muestran en la [[#fig0040|figura 8]] . El caso compresible, no mostrado, tiene un comportamiento similar al caso bidimensional. Nuevamente, se observa el pobre rendimiento del elemento tetraédrico estándar en el campo de los desplazamientos debido al efecto de bloqueo volumétrico. El desplazamiento predicho por este elemento, incluso utilizando mallas tupidas, está muy por debajo del valor correcto. Nótese otra vez que el error en el periodo en la formulación irreducible es notoriamente mayor que con el obtenido de la presente formulación para una misma malla de elementos finitos. Las curvas de convergencia para el caso tridimensional se muestran en la [[#fig0045|figura 9]] . Sobresale la capacidad de MEX-FEM para predecir el valor correcto del desplazamiento vertical (1.998 mm) y la presión (1.251 MPa).

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr8.jpg|center|px|Membrana de Cook 3D. Formulación Irreducible vs MEX-FEM con ν=0.499.]]


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

Figura 8.

Membrana de Cook 3''D'' . Formulación Irreducible vs MEX-FEM con ''ν''  = 0.499.                  

</span>
|}

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr9.jpg|center|px|Membrana de Cook 3D. Curva de convergencia.]]


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

Figura 9.

Membrana de Cook 3''D'' . Curva de convergencia.                  

</span>
|}

Finalmente, la [[#fig0050|figura 10]]  muestra el campo de las presiones de las formulaciones mencionadas en la Malla ''D''  para un valor de ''ν''  = 0.499. Adviértase como MEX-FEM es capaz de dar una solución estabilizada del campo de las presiones en un problema tridimensional.

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

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


[[Image:draft_Content_520834022-1-s2.0-S0213131515000619-gr10.jpg|center|503px|Membrana de Cook 3D. Campo de presiones. Malla D.]]


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

Figura 10.

Membrana de Cook 3''D'' . Campo de presiones. Malla D.                  

</span>
|}

==4. Conclusión==

Ese trabajo propone una formulación mixta explícita de elementos finitos en desplazamientos y deformaciones estabilizadas con sub-escalas, aplicable a elementos con igual orden de interpolación, y en particular, con interpolaciones lineales. Esta formulación se ha desarrollado en el marco del método de las sub-escalas ortogonales, eludiendo de forma efectiva la condición de estabilidad LBB sobre las formulaciones mixtas. En comparación con la formulación irreducible y la formulación mixta sin sub-escala de desplazamientos, la formulación propuesta es más precisa y robusta, no requiere la resolución de un sistema acoplado de ecuaciones y es capaz de abordar problemas quasi-estáticos y dinámicos en situaciones quasi-incompresibles. Por otra parte, el esquema de resolución es análogo al de la formulación irreducible y es muy fácil de integrar la formulación propuesta en un código explícito.

==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 y a los profesores M.W Yuan y Chen Pu por su generoso y valiosísimo apoyo. Riccardo Rossi agradece a los fondos donados por ''The Seventh Framework Programme (FP7/2007-2013)''  del ERC bajo el acuerdo ''n''<sup>''o''</sup>  611636 (NUMEXAS) para el desarrollo de este proyecto. Agradecimientos al programa “Excellence Programme for Knowledge Generation by MINECO” 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. Babuska; Error-bounds for finite element method; Numerische Mathematik, 16 (4) (1971), pp. 322–333</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] 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='bib0015'></span>
[[#bib0015|[3]]] 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='bib0020'></span>
[[#bib0020|[4]]] 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='bib0025'></span>
[[#bib0025|[5]]] M. Chiumenti, M. Cervera, R. Codina; A mixed three-field FE formulation for stress accurate analysis including the incompressible limit; Computer Methods in Applied Mechanics and Engineering, 283 (0) (2015), pp. 1095–1116</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] M. Chiumenti, Q. Valverde, C. Agelet de Saracibar, M. Cervera; A stabilized formulation for incompressible elasticity using linear displacement and pressure interpolations; Computer Methods in Applied Mechanics and Engineering, 191 (46) (2002), pp. 5253–5264</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] R.W. Clough, J. Penzien; Dynamics of structures; Earthquake Engineering & Structural Dynamics (1976)</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] 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='bib0045'></span>
[[#bib0045|[9]]] 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='bib0050'></span>
[[#bib0050|[10]]] R. Codina; Analysis of a stabilized finite element approximation of the oseen equations using orthogonal subscales; Appl. Numer. Math., 58 (3) (2008), pp. 264–283</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] R. Codina, J. Blasco; A finite element method for the stokes problem allowing equal velocity-pressure interpolations; Comput. Meth Appl Mech Eng, 143 (1997), pp. 373–391</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] R. Codina, J. Principe, O. Guasch, S. Badia; Time dependent subscales in the stabilized finite element approximation of incompressible flow problems; Computer Methods in Applied Mechanics and Engineering, 196 (21–24) (2007), pp. 2413–2430</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] 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='bib0070'></span>
[[#bib0070|[14]]] 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='bib0075'></span>
[[#bib0075|[15]]] F. Brezzi, M. Fortin, D. Marini; Mixed finite element methods; Springer (1991)</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] GiD. The personal pre and post processor. http://gid.cimne.upc.es. 2009.</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] A.J. Gil, C.H. Lee, J. Bonet, M. Aguirre; Stabilised petrov-galerkin formulation for linear tetrahedral elements in compressible, nearly incompressible and truly incompressible fast dynamics; Computer Methods in Applied Mechanics and Engineering, 276 (0) (2014), pp. 659–690</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] T.J.R. Hughes; Equivalence of finite elements for nearly incompressible elasticity; J. Appl. Mech., 44 (1977)</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] T.J.R. Hughes; Generalization of selective integration procedures to anisotropic and nonlinear media; International Journal for Numerical Methods in Engineering, 15 (9) (1980), pp. 1413–1418</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]]] T.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]]] T.J.R. Hughes, G.R. Feijóo, L. Mazzei, J.-B. Quincy; The variational multiscale method-a paradigm for computational mechanics; Computer Methods in Applied Mechanics and Engineering, 166 (1–2) (1998), pp. 3–24 Advances in Stabilized Methods 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]]] D.S. Malkus, T.J.R. Hughes; Mixed finite element methods - reduced and selective integration techniques: A unification of concepts; Computer Methods in Applied Mechanics and Engineering, 15 (1) (1978), pp. 63–81</li>
<li><span id='bib0125'></span>
[[#bib0125|[25]]] E. Oñate, A. Valls, J. García; Fic/fem formulation with matrix stabilizing terms for incompressible flows at low and high reynolds numbers; Computational Mechanics, 38 (4–5) (2006), pp. 440–455</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] 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='bib0135'></span>
[[#bib0135|[27]]] G. Scovazzi; Lagrangian shock hydrodynamics on tetrahedral meshes: A stable and accurate variational multiscale approach; Journal of Computational Physics, 231 (24) (2012), pp. 8029–8069</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] J.C. Simo, R.L. Taylor, K.S. Pister; Variational and projection methods for the volume constraint in finite deformation elasto-plasticity; Computer Methods in Applied Mechanics and Engineering, 51 (1–3) (1985), pp. 177–208</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] R.L. Taylor; A mixed-enhanced formulation tetrahedral finite elements; International Journal for Numerical Methods in Engineering, 47 (1-3) (2000), pp. 205–227</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] 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ía, 21 (4) (2005)</li>
</ol>