Triangular element for two-dimensional elasticity with z stress and z strain constant and different from zero.



<big>Elemento triangular para elasticidad bidimensional con &#x03c3;<sub>z</sub>''' ''' y &#x03b5;<sub>z</sub>''' '''constantes y diferentes de cero.</big>

Francisco Zárate <sup>a,b</sup> , Núria Zárate <sup>c </sup>

<sup>a</sup><span style="text-align: center; font-size: 75%;">Centre Internacional de Mètodes Numèrics en Enginyeria (CIMNE) [http://www.cimne.com www.cimne.com]</span>

<sup>b</sup><span style="text-align: center; font-size: 75%;">Universitat Politècnica de Catalunya (UPC),</span>

<span style="text-align: center; font-size: 75%;">Gran Capitan s/n, Campus Nord, 08034 Barcelona, España.<br/></span><sup>c</sup><span style="text-align: center; font-size: 75%;">Universitat de Barcelona (UB)<br/>Gran Via de les Corts Catalanes, 585, 08007 Barcelona, España.</span>

[mailto:zarate@cimne.upc.edu zarate@cimne.upc.edu]

==Resumen==

Este trabajo describe como obtener un elemento bidimensional de tres nodos a partir de un elemento sólido prismático de seis nodos, suponiendo que la tensión normal a las caras triangulares es constante y diferente de cero.

Se presenta también un ejemplo numérico en donde se comparan los resultados obtenidos con el elemento desarrollado y con elementos tridimensionales, poniendo de manifiesto el buen comportamiento del elemento desarrollado.

'''Palabras clave.''' Método de elementos finitos, Tensión y deformación plana, Mecánica de sólidos.

=1.- Introducción=

La gran utilidad de las simplificaciones de la elasticidad tridimensional en problemas de dos dimensiones queda fuera de toda duda pues la reducción en el número de variables permite resolver este tipo de problemas en un tiempo más que razonable, en comparación con el mismo problema resuelto mediante la teoría tridimensional.

Es la elasticidad plana quien estudia la solución particular del problema elástico general en el conjunto de aplicaciones técnicas donde el estado de tensión-deformación es reducible a problemas planos o bidimensionales partiendo de que el dominio de estudio sea un cuerpo prismático. Concretamente se definen dos estados de elasticidad bidimensional: Tensión plana y Deformación plana [1].

Un sólido prismático con una de sus dimensiones mayor que la otras dos es susceptible de ser simplificado en un estado de deformación plana. Para ello se descompone en rebanadas idénticas y se estudia sobre una de estas rebanadas la distribución de deformaciones como un problema bidimensional usando únicamente dos coordenadas <math display="inline">\left( x,y\right)</math>  para definir la posición de cada punto sobre la rebanada. Una condición esencial es que todas las cargas aplicadas estén contenidas en el plano de cada rebanada, y dichas cargas sean las mismas para todas las rebanadas del sólido. De esta manera se genera la hipótesis que define nula la deformación normal a cada una de las rebanadas <math display="inline">\left( {\epsilon }_{z}=\right. </math><math>\left. 0\right)</math> .

Por otra parte, el estado de Tensión plana se puede aplicar a un sólido prismático en el cual una de sus dimensiones es menor que las otras dos, y las cargas aplicadas se encuentran contenidas en la sección media cuya normal es paralela al lado más corto. En este caso, la tensión en dirección a dicha normal se considera nula <math display="inline">\left( {\sigma }_{z}=\right. </math><math>\left. 0\right)</math> .

Sin embargo, estas dos hipótesis no llegan a cubrir la totalidad de los problemas que se pueden simplificar a dos dimensiones. Por ejemplo, el estudio del estado tensional del terreno de un pozo vertical a cierta profundidad como se muestra en la Figura 1, en la cual tanto la deformación como la tensión en dirección vertical son constantes pero diferentes de cero. En este caso es posible definir una rebanada normal a la dirección de perforación del pozo, pero debido a la existencia de una presión normal en la dirección de la rebanada elegida, no es posible aplicar las hipótesis de tensión plana, ni las de deformación plana.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Zarate_347378952-image1-c.png|456px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 1. Aplicación bidimensional donde </span> <math display="inline">{\sigma }_{z}\not =0</math></div>

En el ejemplo mostrado en la Figura 1 se puede considerar que la tensión normal es diferente de cero, pero constante en todo el dominio bidimensional.

Para poder reducir este problema a un estado bidimensional se parte de un elemento finito prismático asociado a un estado de elasticidad tridimensional en donde se aplican las hipótesis de tensión y deformación normal constantes y diferentes de cero.

=2.- Formulación del elemento prismático. =

==2.1-  Hipótesis iniciales==

Como se ha comentado, la formulación del elemento parte del estado tensional tridimensional, en donde se consideran las siguientes hipótesis:

:1. El dominio de análisis corresponde a una laja del dominio de espesor h.

:2. La deformación <math display="inline">{\epsilon }_{z}\,</math> es constante en cualquier punto y diferente de cero.

:3. La tensión <math display="inline">{\sigma }_{z}\,</math> es constante en cualquier punto y diferente de cero.

Como corolario de estas hipótesis se sigue que:

4. Los movimientos en el plano ''xy'' son constantes a lo largo de la posición ''z''.

Estas hipótesis se aplican sobre un elemento prismático tridimensional de seis nodos como el mostrado en la Figura 2.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;"> [[Image:Draft_Zarate_347378952-image2.png|306px]] </span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 2. Elemento prismático base.</span></div>

De manera que siguiendo las hipótesis 2 y 3 se puede definir el estado tensional y de deformaciones en el sólido considerado:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\sigma }_{x}\not =0\, ;\, {\sigma }_{y}\not =0\, ;\, {\sigma }_{z}=cte\, ;\, {\tau }_{xy}\not =0\, \, ;</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\, \, {\tau }_{xz}=0\, ;\, {\tau }_{yz}=0\, \,</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(1)</span>
|}


{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
|  style="text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\epsilon }_{x}\not =0\, ;\, {\epsilon }_{y}\not =0\, ;\, {\epsilon }_{z}=</math><math>cte\, ;\, {\gamma }_{xy}\not =0\, \, ;</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\, \, {\gamma }_{xz}=0\, ;\, {\gamma }_{yz}=0\, \,</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(2)</span>
|}


por lo que describiendo las ecuaciones (1) y (2) en forma vectorial se obtienen las ecuaciones (3) y (4) respectivamente.

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\, \, {\boldsymbol{\sigma }}^{T}=\left[ \begin{matrix}{\sigma }_{x}&{\sigma }_{y}&{\sigma }_{z}&{\tau }_{xy}\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(3)</span>
|}


{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\, {\epsilon }^{T}=\left[ \begin{matrix}{\epsilon }_{x}&{\epsilon }_{y}&{\epsilon }_{z}&{\gamma }_{xy}\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(4)</span>
|}


Y de acuerdo con la teoría de la elasticidad tridimensional, la relación tensión deformación entre (3) y (4) se expresa por:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left[ \begin{matrix}\begin{matrix}{\sigma }_{x}\\{\sigma }_{y}\\{\sigma }_{z}\end{matrix}\\\begin{matrix}{\tau }_{xy}\end{matrix}\end{matrix}\right] =</math><math>\left[ \begin{matrix}{d}_{11}&{d}_{12}&{d}_{12}&0\\{d}_{12}&{d}_{11}&{d}_{12}&0\\{d}_{12}&{d}_{12}&{d}_{11}&0\\0&0&0&{d}_{44}\end{matrix}\right] \left[ \begin{matrix}\begin{matrix}{\epsilon }_{x}\\{\epsilon }_{y}\\{\epsilon }_{z}\end{matrix}\\\begin{matrix}{\gamma }_{xy}\end{matrix}\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(5)</span>
|}


siendo

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{d}_{11}=\frac{E\left( 1-\upsilon \right) }{\left( 1+\upsilon \right) \left( 1-2\upsilon \right) }\, ;\, \, {d}_{12}=</math><math>\frac{E\upsilon }{\left( 1+\upsilon \right) \left( 1-2\upsilon \right) }\quad ;\, {d}_{44}=</math><math>\frac{E}{2\left( 1+\upsilon \right) }\, \,</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(6)</span>
|}

==2.2- Discretización ==

De la hipótesis 4 es posible establecer las siguientes relaciones entre los desplazamientos nodales del elemento:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{u}_{1}={u}_{{1}^{'}}\, ;\, {v}_{1}={v}_{{1}^{'}}\, \, ;</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{u}_{2}={u}_{{2}^{'}}\, ;\, {v}_{2}={v}_{{2}^{'}}\, ;</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{u}_{3}={u}_{{3}^{'}}\, ;\, {v}_{3}={v}_{{3}^{'}}</math><span style="text-align: center; font-size: 75%;"> ;</span>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(7)</span>
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{w}_{1}={w}_{2}={w}_{3}</math><span style="text-align: center; font-size: 75%;"> ;</span>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{w}_{{1}^{'}}={w}_{{2}^{'}}={w}_{{3}^{'}}\,</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|}


Ahora bien, se propone usar un campo bilineal ( <math display="inline">\varnothing =</math><math>a+bx+cy</math>) en el plano <math display="inline">x,\, y\,</math>  para definir el movimiento de los puntos en dicho plano y, como se ha mencionado, dicho plano es constante a lo largo del eje <math display="inline">\, z</math>; por lo que a largo de dicho eje se define un campo lineal e independiente del anterior  ( <math display="inline">\varnothing =</math><math>d+ez</math>). Siendo las funciones de forma que describen el campo bilineal las mostradas en las ecuaciones (8) y las descritas para el eje z las presentadas en las ecuaciones (9).

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{N}_{1}=\frac{1}{2A}\left( {a}_{1}+{b}_{1}x+{c}_{1}y\right)</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{N}_{2}=\frac{1}{2A}\left( {a}_{2}+{b}_{2}x+{c}_{2}y\right)</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(8)</span>
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{N}_{3}=\frac{1}{2A}\left( {a}_{3}+{b}_{3}x+{c}_{3}y\right)</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{N}_{A}=\frac{1}{h}\left( h-z\right)</math> <span style="text-align: center; font-size: 75%;">  ;  </span> <math display="inline">{N}_{B}=</math><math>\frac{1}{h}\left( z\right)</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(9)</span>
|}


Donde <math display="inline">A</math> es el área de la base del elemento prismático de la Figura 2 y <math display="inline">h\,</math>  su altura. Las variables <math display="inline">{a}_{i},\, {b}_{i}</math> y <math display="inline">{c}_{i\, }</math>vienen dadas por las expresiones siguientes, donde los índices <math display="inline">i,\, j,\, k</math> son cíclicos.

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{a}_{i}={x}_{j}{y}_{k}-{x}_{k}{y}_{j}\,</math> 
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{b}_{i}={y}_{j}-{y}_{k}</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(10)</span>
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{c}_{i}={x}_{k}-{x}_{j}\, \, ;\, \, i,j,k=1,2,3</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|}


Si bien, podría bastar un campo constante a lo largo del eje z en lugar de uno lineal como se propone, se ha optado por utilizar un campo lineal a fin de mantener lo más homogénea la descripción de las variables nodales.

Las deformaciones definidas en la ecuación (4) quedarán definidas como:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\epsilon }_{x}=\, \frac{\partial {N}_{1}}{\partial x}{u}_{1}+\frac{\partial {N}_{2}}{\partial x}{u}_{2}+</math><math>\frac{\partial {N}_{3}}{\partial x}{u}_{3}\, \, ;\, {\epsilon }_{y}=\, \frac{\partial {N}_{1}}{\partial y}{v}_{1}+</math><math>\frac{\partial {N}_{2}}{\partial y}{v}_{2}+\frac{\partial {N}_{3}}{\partial y}{v}_{3}</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| 
|}
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\epsilon }_{z}=\, \frac{\partial {N}_{A}}{\partial z}{w}_{1}+\frac{\partial {N}_{B}}{\partial z}{w}_{{1}^{'}}</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(11)</span>
|-
|  style="text-align: center;text-align: right;white-space: nowrap;"|
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\gamma }_{xy}=\, \frac{\partial {N}_{1}}{\partial x}{v}_{1}+\frac{\partial {N}_{2}}{\partial x}{v}_{2}+</math><math>\frac{\partial {N}_{3}}{\partial x}{v}_{3}+\frac{\partial {N}_{1}}{\partial y}{u}_{1}+</math><math>\frac{\partial {N}_{2}}{\partial y}{u}_{2}+\frac{\partial {N}_{3}}{\partial y}{u}_{3}</math>
|}
|  style="text-align: right;text-align: right;white-space: nowrap;"|
|}


Y la descripción del vector de desplazamientos viene dada por:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\boldsymbol{a}}^{T}=\left[ \begin{matrix}{u}_{1}&{v}_{1}&{u}_{2}&{v}_{2}&{u}_{3}&{v}_{3}&{w}_{1}&{w}_{{1}^{'}}\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(12)</span>
|}


Por lo que, haciendo uso de las funciones de forma descritas, la matriz '''B''' de deformaciones queda definida como:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\boldsymbol{B}=\left[ \begin{matrix}\frac{{b}_{1}}{2A}&0&\frac{{b}_{2}}{2A}&0&\frac{{b}_{3}}{2A}&0&0&0\\0&\frac{{c}_{1}}{2A}&0&\frac{{c}_{2}}{2A}&0&\frac{{c}_{3}}{2A}&0&0\\0&0&0&0&0&0&-\frac{1}{h}&\frac{1}{h}\\\frac{{c}_{1}}{2A}&\frac{{b}_{1}}{2A}&\frac{{c}_{2}}{2A}&\frac{{b}_{2}}{2A}&\frac{{c}_{3}}{2A}&\frac{{b}_{3}}{2A}&0&0\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(13)</span>
|}


De aquí es muy sencillo obtener la matriz de rigidez del elemento aplicando el principio de los trabajos virtuales:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\boldsymbol{K}}^{(e)}=\int_{{A}^{(e)}}^{}{\boldsymbol{B}}^{T}\boldsymbol{DB}hd{A}^{(e)}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(14)</span>
|}

==2.3- Simplificación de variables ==

La matriz de rigidez mostrada en la ecuación (14) es de <math display="inline">8\times 8</math> y responde a las variables nodales definidas en la ecuación (12) por lo que es necesario eliminar las variables <math display="inline">{w}_{1}</math><span style="text-align: center; font-size: 75%;"> y </span> <math display="inline">{w}_{{1}^{'}}.</math>

De la relación tensión deformación para la elasticidad 3D, mostrada en la ecuación (5) y dado que se considera que <math display="inline">{\sigma }_{z}=</math><math>cte</math> en todo el dominio, se puede escribir:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\sigma }_{z}=\, {d}_{12}\, {\epsilon }_{x}+\, {d}_{12}\, {\epsilon }_{y}+\, {d}_{11}\, {\epsilon }_{z}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(15)</span>
|}


de donde se desprende que:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\epsilon }_{z}=\, \frac{1}{{d}_{11}}\left( {{\sigma }_{z}-d}_{12}\, {\epsilon }_{x}-\right. </math><math>\left. \, {d}_{12}\, {\epsilon }_{y}\, \right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(16)</span>
|}


y por lo tanto, la elongación del segmento <math display="inline">\overline{{\, w}_{1}\, {w}_{{1}^{'}}}{=d}_{z}</math>que se define como

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{d}_{z}=h{\epsilon }_{z}=\, \frac{h}{{d}_{11}}\left( {{\sigma }_{z}-d}_{12}\, {\epsilon }_{x}-\right. </math><math>\left. \, {d}_{12}\, {\epsilon }_{y}\, \right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(17)</span>
|}


Por lo que si el desplazamiento  <math display="inline">{w}_{1}=0</math> entonces

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{w}_{{1}^{'}}=\frac{h}{{d}_{11}}\left( {{\sigma }_{z}-d}_{12}\, {\epsilon }_{x}-\right. </math><math>\left. \, {d}_{12}\, {\epsilon }_{y}\, \right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(18)</span>
|}


Obsérvese ahora que las dos variables asociadas al espesor del elemento han sido eliminadas pues <math display="inline">{w}_{1}=</math><math>0</math> y <math display="inline">{w}_{{1}^{'}}</math> es ahora función de <math display="inline">{\epsilon }_{x},\, {\epsilon }_{y}\,</math> y <math display="inline">{\sigma }_{z}</math> por lo que se define un nuevo vector de variables nodales <math display="inline">{\boldsymbol{a}}_{1}</math> como:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\boldsymbol{a}}_{1}^{T}=\left[ \begin{matrix}{u}_{1}&{v}_{1}&{u}_{2}&{v}_{2}&{u}_{3}&{v}_{3}&{\sigma }_{z}\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(19)</span>
|}


De manera que la relación entre <math display="inline">{\boldsymbol{a}}_{1}</math> y <math display="inline">\boldsymbol{a}</math> queda definida por:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\boldsymbol{a}=\boldsymbol{V}{\boldsymbol{a}}_{\mathit{\boldsymbol{1}}}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(20)</span>
|}


en donde

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\boldsymbol{V}=\left[ \begin{matrix}1&0&0&0&0&0&0\\0&1&0&0&0&0&0\\0&0&1&0&0&0&0\\0&0&0&1&0&0&0\\0&0&0&0&1&0&0\\0&0&0&0&0&1&0\\0&0&0&0&0&0&0\\\frac{-{d}_{12}{b}_{1}h}{{2Ad}_{11}}&\frac{-{d}_{12}{c}_{1}h}{2A{d}_{11}}&\frac{-{d}_{12}{b}_{2}h}{2A{d}_{11}}&\frac{-{d}_{12}{c}_{2}h}{{2Ad}_{11}}&\frac{-{d}_{12}{b}_{3}h}{2A{d}_{11}}&\frac{-{d}_{12}{c}_{2}h}{2A{d}_{11}}&\frac{h}{{d}_{11}}\end{matrix}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(21)</span>
|}


Al utilizar la definición de <math display="inline">{\boldsymbol{a}}_{1}</math>como variables nodales elementales en la ecuación de equilibrio se tiene:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\boldsymbol{K}}_{\mathit{\boldsymbol{v}}}^{\mathit{\boldsymbol{(e)}}}{\boldsymbol{a}}_{\mathit{\boldsymbol{1}}}^{\mathit{\boldsymbol{(e)}}}=</math><math>{\boldsymbol{f}}^{\mathit{\boldsymbol{(e)}}}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(22)</span>
|}


en donde

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\boldsymbol{K}}_{\mathit{\boldsymbol{v}}}^{\mathit{\boldsymbol{(e)}}}={\boldsymbol{K}}_{}^{\mathit{\boldsymbol{(e)}}}{\boldsymbol{V}}^{\mathit{\boldsymbol{(e)}}}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span style="text-align: center; font-size: 75%;">(23)</span>
|}


La matriz de rigidez <math display="inline">{\boldsymbol{K}}_{\mathit{\boldsymbol{v}}}^{(\mathit{\boldsymbol{e}})}</math> que se obtiene es de <math display="inline">8\times 7</math>. Sin embargo, la séptima ecuación del sistema mostrado en (22) puede ser eliminar directamente ya que resulta ser una ecuación linealmente dependiente de las demás y está asociada al desplazamiento <math display="inline">{w}_{1}\,</math>  cuyo valor es conocido e igual a 0, quedando un sistema de <math display="inline">7\times 7</math>.

Finalmente, dado que <math display="inline">\, {\sigma }_{z}</math> no es una incógnita sino un valor conocido y constante en todo el domino, es posible eliminar la ecuación asociada, actualizando el vector <math display="inline">{\boldsymbol{f}}^{(\mathit{\boldsymbol{e}})}</math> con el producto de <math display="inline">{\, K}_{i,7\, }\ast \, {\sigma }_{z}</math> obteniendo una matriz de 6 x 6 dependiente exclusivamente de los desplazamientos en el plano <math display="inline">x,\, y\, \,</math> del dominio por analizar.

El ensamblaje y solución del problema se sigue de la manera estándar a cualquier problema bidimensional, como se muestra en el siguiente apartado.

Esta formulación se ha implementado en el código FEM2DEM desarrollado en CIMNE y siguiendo la formulación propuesta en [2-4].

=3.- Ejemplos. =

En este apartado se muestran dos ejemplos, el primero corresponde a la verificación de los desplazamientos en un dominio rectangular debido a la existencia de un esfuerzo <math display="inline">{\sigma }_{z}</math>. Los resultados se comparan con los obtenidos mediante un programa de solidos tridimensionales a fin de verificar que los resultados obtenidos sean iguales.

El segundo ejemplo corresponde a una sección de un pozo de perforación en donde además de considerar la profundidad se evalúa el efecto de la presión interna del ducto de perforación.

==3.1- Ejemplo de validación. ==


[[Image:Draft_Zarate_347378952-image3-c.png|center|498px]]
Este ejemplo corresponde a un cubo de dimensiones <math display="inline">1.2\times 1.2\times 1.2\, m.</math>  el cual está sujeto a una presión axial de <math display="inline">1\times {10}^{7\, }\, N/{m}^{2}</math> en compresión y donde las propiedades del material son: <math display="inline">E=</math><math>2.1\times {10}^{11}\, N/{m}^{2}</math>, <math display="inline">\nu =0.20\,</math> . En la Figura 3 se muestran las condiciones de carga (Fig. 3 a) y de apoyo (Fig. 3 b) el cual restringe, sobre ambas superficies de carga, el movimiento del punto medio de las aristas en dirección normal. Por otra parte, no es necesario restringir el desplazamiento vertical dado que las cargas de compresión son auto equilibradas. Se ha procurado que ambas mallas tengan una cierta similitud en cuanto al tamaño del elemento para evitar diferencias inducidas por la calidad de la malla.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 3. Ejemplo de validación. a) aplicación de cargas. b) restricción de movimientos.</span></div>

En la Figura 4 se muestran los resultados del análisis tridimensional los cuales son usados de referencia. Se observa que tanto el desplazamiento en <math display="inline">x\,</math>  como en <math display="inline">\, y\, \,</math> son idénticos y con valores totalmente simétricos. Siendo el desplazamiento máximo de <math display="inline">1.1406\times {10}^{-5}\, m</math>.

La Figura 5 muestra los resultados de los mismos desplazamientos sobre un dominio bidimensional utilizando la formulación propuesta. En este caso, también los desplazamientos en <math display="inline">x</math> e <math display="inline">y</math> son idénticos y simétricos con un valor máximo de <math display="inline">1.1429\times {10}^{-5}\, m</math>.

Comparando ambos análisis se puede concluir que la formulacion presentada es congruente y correcta ya que el error entre ambas formulaciones es prácticamente nulo ( <math display="inline">0.20\%</math> ).

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Zarate_347378952-image4-c.png|center|516px]]

[[Image:Draft_Zarate_347378952-image5-c.png|center|528px]]
<span style="text-align: center; font-size: 75%;">Figura 4.  Ejemplo de validación. Resultados 3D sobre la superficie de carga de los desplazamientos </span> <math display="inline">x</math><span style="text-align: center; font-size: 75%;"> e </span> <math display="inline">y</math><span style="text-align: center; font-size: 75%;">.</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 5.  Ejemplo de validación. Resultados 2D de los desplazamientos </span> <math display="inline">x</math><span style="text-align: center; font-size: 75%;"> e </span> <math display="inline">y</math></div>

==3.2- Análisis triaxial. ==

En este ejemplo se reproduce un estado triaxial sobre un cilindro horadado y de material homogéneo por el cual circula un fluido a una presión similar a la exterior. Axialmente el cilindro se encuentra comprimido a <math display="inline">1.0\times {10}^{7}\, Pa</math> mientras que se aplica una presión sobre la superficie exterior del cilindro de <math display="inline">5.0\times {10}^{6}\, Pa</math> al tiempo que la presión interna es de <math display="inline">1.0\times {10}^{7}\, Pa</math>. Este ejemplo reproduce una sección de un pozo sin camisa por donde se inyecta un fluido a presión. El análisis clásico es tridimensional, sin embargo, aplicando la formulación presentada es posible atacar el problema mediante un análisis 2D.

<span id='_GoBack'></span>La probeta se encuentra formada por un disco de diámetro exterior de <math display="inline">0.1016\, m</math>. y de diámetro interior de <math display="inline">0.0254\, m</math> como se muestra en la Figura 6 a). El material a usar tiene un Módulo de elasticidad de <math display="inline">3.4\times {10}^{9}\, Pa</math> y coeficiente de Poisson de <math display="inline">0.20</math>.

La malla utilizada se muestra en la Figura 6 b) y está formada por 352 nodos y 616 elementos triangulares de tres nodos.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Zarate_347378952-image6-c.png|474px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 6.  Análisis triaxial. a) Geometría. b) Malla utilizada.</span></div>

Los resultados obtenidos se muestran en la Figura 7 en donde se puede observar los desplazamientos tanto en áreas de igual color (a) como de forma vectorial (b).

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Zarate_347378952-image7-c.png|522px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 7.  Análisis triaxial. Resultados 2D de los desplazamientos </span> <math display="inline">x</math><span style="text-align: center; font-size: 75%;"> e </span> <math display="inline">y</math></div>

La Figura 8 muestra la distribución de todos los esfuerzos involucrados en la formulación y quizás este sea el aspecto más interesante pues se puede describir un comportamiento tridimensional de forma adecuada, sobre todo, en el caso ser necesario un análisis plástico o de daño. La Figura 7 a) y 7 b) representan los esfuerzos <math display="inline">\, {\sigma }_{x}</math> y <math display="inline">\, {\sigma }_{y}</math> respectivamente. La Figura 7 c) corresponde al esfuerzo <math display="inline">\, {\sigma }_{z}\,</math> y finalmente, la figura 7 d) el esfuerzo de corte <math display="inline">\, {\sigma \gamma }_{xy}</math>.

==Conclusiones==

El elemento presentado tiene una buena respuesta y es capaz de evaluar adecuadamente el efecto de tener un esfuerzo normal diferente de cero. Una de las ventajas asociadas a esta formulación es la de poder modelar mediante un problema bidimensional el daño o plasticidad en un estado triaxial sin necesidad de tener que resolver un problema 3D, con el consabido error de tiempo y memoria.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Zarate_347378952-image8-c.png|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figura 8.  Análisis triaxial. a) esfuerzos </span> <math display="inline">\, {\sigma }_{x}</math><span style="text-align: center; font-size: 75%;"> b) esfuerzos </span> <math display="inline">\, {\sigma }_{y}</math><span style="text-align: center; font-size: 75%;"> c) esfuerzos </span> <math display="inline">\, {\sigma }_{z}</math><span style="text-align: center; font-size: 75%;"> d) esfuerzos </span> <math display="inline">\, {\gamma }_{xy}</math></div>

==Referencias==

[1] Oñate, E. ''Cálculo de Estructuras por el Método de los Elementos Finitos. Análisis Estático Lineal''. Vol 1. (2016) ISBN: 978-84-945689-7-8 Editorial: CIMNE Barcelona. España

<span style="text-align: center; font-size: 75%;">[2]</span> <span style="text-align: center; font-size: 75%;">Zárate, F. & Oñate, E. A simple FEM–DEM technique for fracture prediction in materials and structures.  Comp. Part. Mech. (2015) 2: 301. doi.org/10.1007/s40571-015-0067-2</span>

<span style="text-align: center; font-size: 75%;">[3]</span> <span style="text-align: center; font-size: 75%;">Zárate, F., Cornejo, A. & Oñate E., A three-dimensional FEM–DEM technique for predicting the evolution of fracture in geomaterials and concrete. Comp. Part. Mech. (2018) 5: 411. doi.org/10.1007/s40571-017-0178-z</span>

<span style="text-align: center; font-size: 75%;">[4]</span> <span style="text-align: center; font-size: 75%;">Zárate, F. & Oñate E. Predicción de fracturas en estructuras de hormigón combinando los métodos de elementos finitos y de elementos discretos. Hormigón y Acero. (2018) doi.org/10.1016/j.hya.2018.05.002 </span>