Matias RIMNI

==Resumen==

En este trabajo se describe la formulación corrotacional de un elemento de viga que utiliza distintas teorías, Euler-Bernoulli, Timoshenko y Euler-Bernoulli incluyendo el acoplamiento entre los efectos de los esfuerzos de flexión y axil lo que genera coeficientes no lineales en la matriz de rigidez elástica de dicho elemento.  La cinemática corrotacional se basa en la separación del movimiento de un sólido  en una parte deformacional y otra en movimiento de cuerpo rígido.  La parte deformacional del movimiento es  representada por tres grados de libertad: una traslación debido al esfuerzo axil y por las rotaciones nodales debidas a la flexión del elemento. A nivel local se obtienen los esfuerzos internos  mediante el princípio de los trabajos virtuales una vez definidas las funciones de interporlación  que describen el movimiento deformacional para cada teoría de viga  utilizada.  Se describe de forma detallada la obtención de las matrices de rigideces  elástica para los tres tipos de elementos de viga.  A través de algunos ejemplos numéricos se muestra la habilidad del elemento para obtener el comportamiento no lineal de vigas, pórticos y arcos.

'''Palabras claves''': Elemento de viga Euler-Bernoulli. Elemento de viga Timoshenko,  Elemento de viga Euler-Bernoulli no lineal, Formulación corrotacional.

==Abstract==

The present work uses a co-rotating approach for obtain the internal forces and tangent stiffness matrices for three plane beam elements. All of them are based on the same co-rotational approach, and differ by the strain definition used in the local co-rotational coordinate system. Based on the Bernoulli assumption, the first two elements use a linear and a shallow arch strain definition, respectively. The third element is based on the Timoshenko assumption with linear interpolations for the displacements. Some examples are presented and the numerical results demonstrate that the beam elements here presented are able to perform the nonlinear analysis of plane frames and 2D arches and to handle large rotations too.

==1 Introducción==

En las últimas dos décadas se ha observado un creciente interés en la comunidad de mecánica computacional referente a la aplicación de la formulación corrotacional en el análisis no lineal de estructuras utilizadas en la industria aeronáutica, aeroespacial y petrolera. Felippa y Haugen [1] describen en detalle una reseña del desarrollo de la formulación corrotacional y de sus aplicaciones en mecánica computacional. La idea principal de la formulación corrotacional es descomponer el movimiento de un sólido en la suma de un movimiento de cuerpo rígido y un movimiento deformacional. Esta formulación es un caso particular de la cinemática Lagrangiana para modelar los efectos de la no linealidad geométrica en el análisis estructural, haciendo uso del método de los elementos finitos. En esta formulación se asume que tanto las traslaciones como las rotaciones de cuerpo rígido podrán ser grandes, pero las deformaciones han de ser pequeñas y por ello existe la posibilidad del uso de elementos finitos lineales en problemas que involucran no linealidad geométrica, siendo esta la principal motivación del uso de la formulación corrotacional. En el panorama actual se puede citar algunos trabajos que han aplicado la formulación corrotacional utilizando otros métodos de discretización, por ejemplo, en Yaw et al. [2] se muestran aplicaciones de la formulación corrotacional y los métodos sin malla en el análisis de sólidos 2D incluyendo no linealidad geométrica y física; Xu et al [3] aplican el método de elementos finitos extendidos para el tratamiento de problemas de localización y la formulación corrotacional en elementos de viga de Timoshenko en el análisis no lineal de estructuras reticulares. Aún en el contexto de la aplicación de la formulación corrotacional en la actualidad, en Bisegna y Caselli [4] se describe la formulación corrotacional de un elemento de lámina hiperelástico y su aplicación en el estudio de fibras de tejidos biológicos. Como ejemplo más reciente de aplicación de la formulación corrotacional en el análisis no lineal de estructuras se puede citar el trabajo de Li et al. [5] que describe la formulación corrotacional de um elemento triangular de lámina de 6 nodos y su aplicación en el análisis elastoplastico de cáscaras finas. En este trabajo se aplica la formulación corrotacional junto a diferentes teorías de vigas en el análisis no lineal geométrico de vigas, pórticos y arcos. Para ello, se utilizan las teorías clásicas de viga de Euler-Bernoulli y de Timoshenko. En el tercer elemento de viga se incluye en la teoría de Euler-Bernoulli el efecto del acoplamiento de los esfuerzos axil y de flexión, dicho acoplamiento genera coeficientes no lineales en la matriz de rigidez elástica del elemento, razón por la cual se le denonima elemento de viga de Euler- Bernoulli ''no lineal''. La cinemática corrotacional aquí descrita sigue la formulación propuesta  por Crisfield [6], esta formulación permite obtener la matriz de rigidez tangente en coordenadas globales independientemente del tipo de elemento de viga adoptado. El objetivo en este trabajo es comparar el desempeño de los tres tipos de elementos de viga 2D al retratar el comportamiento no lineal de vigas, pórticos planos y arcos.

==2 Cinemática corrotacional==

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_709676212-corotating_2Dbeam_element.png|386px|Movimiento del elemento de viga 2D.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 1:''' Movimiento del elemento de viga 2D.
|}

La Figura 1 muestra las variables cinemáticas utilizadas en la descripción corrotacional del elemento de viga 2D. Las coordenadas nodales del elemento en la configuración de referencia son <math display="inline">(x_{_1},y_{_1})</math> y <math display="inline">(x_{_2},y_{_2})</math>, respectivamente. En la configuración deformada dichas coordenadas se definem como <math display="inline">{(x'}_{_1}{,y'}_{_1})</math> y <math display="inline">{(x'}_{_2}{,y'}_{_2}).</math> A cada elemento se le asigna el sistema de referencia local <math display="inline">(x_{_l},y_{_l})</math> que lo acompaña durante su movimiento. El sistema de referencia local es utilizado para definir la parte deformacional del movimiento. El vector de desplazamientos globales es definido por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf u}_{_g}=(u_{_1}, v_{_1}, \theta _{_1}, u_{_2}, v_{_2}, \theta _{_2})^{^T} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

La parte deformacional del movimiento es dada por el vector de desplazamientos locales como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf u}_{_l}=(\bar u, \bar \theta _{_1}, \bar \theta _{_2})^{^T} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

donde <math display="inline">\bar u</math> es el desplazamiento relativo entre los nodos 1 y 2 en la dirección del eje local <math display="inline">x_{_l}</math>. <math display="inline">\bar \theta _{_1}</math> y <math display="inline">\bar \theta _{_2}</math> son las rotaciones deformacionales de los nodos 1 y 2, respectivamente. Las componentes de <math display="inline">{\bf u}_{_l}</math> son definidas de la siguiente manera

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar u=l-l_{_0}</math>
|-
| style="text-align: center;" | <math>\bar \theta _{_1}=\theta _{_1}-\alpha</math>
|-
| style="text-align: center;" | <math>\bar \theta _{_2}=\theta _{_2}-\alpha</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

donde <math display="inline">l_{_0}</math> y <math display="inline">l</math> denotan la longitud inicial y actual del elemento. Dichas longitudes son dadas por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>l_{_0}=\sqrt{x_{_{21}}^2+y_{_{21}}^2}</math>
|-
| style="text-align: center;" | <math>l=\sqrt{(x_{_{21}}+u_{_{21}})^2+(y_{_{21}}+v_{_{21}})^2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>x_{_{21}}=x_{_2}-x_{_1}</math>
|-
| style="text-align: center;" | <math>y_{_{21}}=y_{_2}-y_{_1}</math>
|-
| style="text-align: center;" | <math>u_{_{21}}=u_{_2}-u_{_1}</math>
|-
| style="text-align: center;" | <math>v_{_{21}}=v_{_2}-v_{_1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

y <math display="inline">\alpha=\beta-\beta_{_0}</math> denota la rotación de cuerpo rígido que puede ser calculada como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{sen}\alpha={cos}\beta _{_0}{sen}\beta{-}{sen}\beta _{_0}{cos}\beta</math>
|-
| style="text-align: center;" | <math>{cos}\alpha={cos}\beta _{_0}{cos}\beta{+}{sen}\beta _{_0}{sen}\beta</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

donde <math display="inline">\beta _{_0}</math> es la inclinación del elemento en la configuración de referencia y <math display="inline">\beta </math> la inclinación de la línea que une los nodos del elemento en la configuración actual. Estos ángulos pueden ser calculados con el uso de las siguientes expresiones

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{cos}\beta _{_0}=\frac{x_{_{21}}}{l_{_0}}</math>
|-
| style="text-align: center;" | <math>{sen}\beta _{_0}=\frac{y_{_{21}}}{l_{_0}}</math>
|-
| style="text-align: center;" | <math>{cos}\beta =\frac{x_{_{21}}+u_{_{21}}}{l}=\frac{{x'}_{_{21}}}{l}</math>
|-
| style="text-align: center;" | <math>{sen}\beta =\frac{y_{_{21}}+v_{_{21}}}{l}=\frac{{y'}_{_{21}}}{l}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

Usando la ecuación (7) se puede reescribir la ecuación (6) como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}{cos}\alpha =\frac{x_{_{21}}^{\ast }}{l_{_0}l},~~~     {sen}\alpha =\frac{y_{_{21}}^{\ast }}{l_{_0}l}     \\     x_{_{21}}^{\ast }=x_{_{21}}{x'}_{_{21}}+y_{_{21}}{y'}_{_{21}}     \\     y_{_{21}}^{\ast }=x_{_{21}}{y'}_{_{21}}-y_{_{21}}{x'}_{_{21}}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

donde <math display="inline">x_{_{21}}^{\ast }</math> y <math display="inline">y_{_{21}}^{\ast }</math> estan relacionadas con las diferencias entre las traslaciones de cuerpo rígido de los nodos del elemento en las direciones <math display="inline">x</math> e <math display="inline">y</math>, respectivamente. Por otra parte, utilizando la relación <math display="inline">{tan}\frac{\alpha }{2}=\frac{1-{cos}\alpha }{{sin}\alpha }</math>,  la rotación de cuerpo rígido es dada por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha=2^{-1}\Bigl(\frac{1-{cos\alpha }}{{sen}\alpha }\Bigr)=2{tan}^{-1}   \Bigl(\frac{l_{_0}l-x_{_{21}}^{\ast }}{y_{_{21}}^{\ast }}\Bigr) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

Esta expresión es singular si <math display="inline">y_{_{21}}^{\ast }=0</math>, por lo que será asignado el valor de <math display="inline">\alpha=0</math> para <math display="inline">x_{_{21}}^{\ast }=l_{_0}l</math> o <math display="inline">\alpha =\pi </math> para <math display="inline">x_{_{21}}^{\ast }=-l_{_0}l</math>. Para el cálculo de las rotaciones deformacionales se utilizan las ecuaciones (3b) y (3c). Con el objeto de asumir cualquier magnitud en el valor de los ángulos, las expresiones para las rotaciones deformacionales deben ser calculadas usando la función módulo [7], tal que:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\bar \theta _{_1}=\underset{2\pi }{mod}(\bar \theta _{_1}+\pi )-\pi      \\     \bar \theta _{_2}=\underset{2\pi }{mod}(\bar \theta _{_2}+\pi )-\pi    \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

La función módulo asigna un valor arbitrario en el intervalo <math display="inline">[0,2\pi [</math>, y el último término restablece la simetría en relación al valor cero. La omisión de este paso puede llevar a problemas de divergencia cuando el elemento de viga sufra rotaciones <math display="inline">\pm \pi ,\pm 2\pi ,...</math>

===2.1 Relación entre variables locales y globales===

Para obtener la relación entre variables locales y globales se aplica una variación virtual en los desplazamientos locales, tal que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \bar u       =\delta l=     \bigl(-{cos}\beta , -{sen}\beta , 0, {cos}\beta , {sen}\beta , 0\bigr)^{^T}\delta{\bf u}_{_g}</math>
|-
| style="text-align: center;" | <math>\delta \bar \theta _{_1}=\delta \theta _{_1}-\delta \alpha =\delta \theta _{_1}-\delta \beta </math>
|-
| style="text-align: center;" | <math>\delta \bar \theta _{_2}=\delta \theta _{_2}-\delta \alpha =\delta \theta _{_2}-\delta \beta</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

La ecuación (11a) se obtiene diferenciando la ecuación (4b) respecto a los desplazamientos globales. Por otra parte, la variación virtual de <math display="inline">\beta </math> es calculada diferenciando la ecuación (3.7d) respecto a los desplazamientos globales y se escribe como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \beta =\frac{1}{l}\bigl({sen}\beta , -{cos}\beta , 0, -{sen}\beta , {cos}\beta , 0\bigr)^{^T}   \delta{\bf u}_{_g} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

Reagrupando en forma matricial las relaciones entre las variaciones virtuales de las variables locales y globales, se escribe la matriz de transformación <math display="inline">{\bf B}</math> como
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size: 100%;">
<math>\delta{\bf u}_{_l}={\bf B}\delta{\bf u}_{_g}</math></div>

con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf B}=   \begin{bmatrix}-{cos}\beta & -{sen}\beta & 0 & {cos}\beta & {sen}\beta & 0 \\     -\frac{{sen}\beta }{l} & \frac{{cos}\beta }{l} & 1 & \frac{{sen}\beta }{l} & -\frac{{cos}\beta }{l} & 0     \\     -\frac{{sen}\beta }{l} & \frac{{cos}\beta }{l} & 0 & \frac{{sen}\beta }{l} & -\frac{{cos}\beta }{l} & 1   \end{bmatrix}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

Por otro lado, para obtener la relación entre los esfuerzos internos y los esfuerzos globales se impone que el trabajo virtual sea igual entre dichos esfuerzos, asi que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V=\delta{\bf u}_{_g}^{\scriptscriptstyle{T}}{\bf f}_{_g}=\delta{\bf u}_{_l}^{\scriptscriptstyle{T}}{f}_{_l}=     \delta{\bf u}_{_g}^{\scriptscriptstyle{T}}{\bf B}^{^T}{\bf f}_{_l} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

La igualdad descrita en la ecuación (14) es válida para cualquier valor arbitrario de <math display="inline">\delta{\bf u}_{_g}</math>, de modo que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf f}_{_g}={\bf B}^{^T}{\bf f}_{_l} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

donde el vector de esfuerzos internos <math display="inline">{\bf f}_{_l}=(N, M_{_1}, M_{_2})^{^T}</math> depende de la definición del elemento de viga utilizado. En los próximos apartados se obtienen dichos vectores.

===2.2 Matriz de rigidez tangente global===

La relación entre las variaciones virtuales del vector de fuerzas global y del vector de desplazamientos global es descrita por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta{\bf f}_{_g}={\bf K}_{_g}\delta{\bf u}_{_g} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

donde <math display="inline">{\bf K}_{_g}</math> es la matriz de rigidez tangente global que se obtiene diferenciando la ecuación (15) como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta{\bf f}_{_g}={\bf B}^{^T}\delta{\bf f}_{_l}+\delta{\bf B}^{^T}{\bf f}_{_l}=   {\bf B}^{^T}\delta{\bf f}_{_l}+N\delta{\bf b}_{_1}+M_{_1}\delta{\bf b}_{_2}+M_{_2}\delta{\bf b}_{_3} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

donde <math display="inline">{\bf b}_{_1}</math>, <math display="inline">{\bf b}_{_2}</math> y <math display="inline">{\bf b}_{_3}</math> son la primera, segunda y tercera columna de la matriz <math display="inline">{\bf B}^{^T}</math>. Para auxiliar en el desarrollo algebraico son definidos los siguientes vectores

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c} r=&\bigl(-{cos}\beta , -{sen}\beta , 0, {cos}\beta , {sen}\beta , 0\bigr)^{^T}     \\     z=&\bigl({sen}\beta , -{cos}\beta , 0, -{sen}\beta , {cos}\beta , 0\bigr)^{^T}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}
<!--
cuyas variaciones virtuales son

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\delta{\bf r}=& {\bf z} \delta \beta \\     \delta{\bf z}=&-{\bf r}\delta \beta    \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

De modo que, las ecuaciones (11a) y (12) pueden reescribirse de la siguiente forma

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\delta \bar u=&\delta l={\bf r}^{\scriptscriptstyle{T}}\delta{\bf u}_{_g}\\     
\delta \beta =&\frac{1}{l}{\bf z}^{\scriptscriptstyle{T}}\delta{\bf u}_{_g}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

Con esta notación los vectores <math display="inline">{\bf b}_{_1}</math>, <math display="inline">{\bf b}_{_2}</math> y <math display="inline">{\bf b}_{_3}</math> se escriben como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}{\bf b}_{_1}=& {\bf r}\\     {\bf b}_{_2}=&(0, 0, 1, 0, 0, 0)^{^T}-\frac{1}{l}{\bf z}\\     {\bf b}_{_3}=&(0, 0, 0, 0, 0, 1)^{^T}-\frac{1}{l}{\bf z}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

cuyas variaciones virtuales son dadas por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\delta{\bf b}_{_1}=&\delta{\bf r}=\frac{1}{l}{\bf z}{\bf z}^{\scriptscriptstyle{T}}\delta{\bf u}_{_g}     \\     \delta{\bf b}_{_2}=&\delta{\bf b}_{_3}=-\frac{1}{l}\delta{\bf z}+\frac{\delta l}{l^{^2}}{\bf z}=     \frac{1}{l^{^2}}({\bf r}{\bf z}^{\scriptscriptstyle{T}}+{\bf z}{\bf r}^{\scriptscriptstyle{T}})\delta{\bf u}_{_g}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

El primer término de la ecuación (17) es calculado al introducir la matriz de rigidez elástica local <math display="inline">{\bf K}_{_l}</math>, que depende del tipo de elemento de viga adoptado. Esta matriz será calculada en los próximos apartados. Por lo tanto, la relación entre la variación virtual de los esfuerzos internos y los desplazamiento globales es definida por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta{\bf f}_{_l}={\bf K}_{_l}\delta{\bf u}_{_l}={\bf K}_{_l}{\bf B}\delta{\bf u}_{_g} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Finalmente, de las ecuaciones (16), (17), (22a), (22b) y (23), se obtiene la expresión de la matriz de rigidez tangente global como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf K}_{_g}={\bf B}^{^T}{\bf K}_{_l}{\bf B}+{\bf z}{\bf z}^{\scriptscriptstyle{T}}\frac{N}{l}+\frac{1}{l^{^2}}   ({\bf r}{\bf z}^{\scriptscriptstyle{T}}+{\bf z}{\bf r}^{\scriptscriptstyle{T}})(M_{_1}+M_{_2}) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

El primer término de la ecuación (24) representa la matriz de rigidez material. Así mismo, teniendo en cuenta que el esfuerzo cortante se define como <math display="inline">Q=\frac{M_{_1}+M_{_2}}{l}</math>, la matriz de rigidez geométrica global puede ser escrita mediante la siguiente expresión

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf K}_{_g}^{\ast }=\frac{N}{l}{\bf z}{\bf z}^{\scriptscriptstyle{T}}+\frac{Q}{l}({\bf r}{\bf z}^{\scriptscriptstyle{T}}+{\bf z}{\bf r}^{\scriptscriptstyle{T}}) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

Por otra parte, la matriz de rigidez elastica <math display="inline">{\bf K}_{_l}</math> es obtenida por medio de la deferenciación de los esfuerzos internos <math display="inline">{\bf f}_{_l}</math> respecto a los desplazamientos locales <math display="inline">{\bf u}_{_l}</math>, tal que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf K}_{_l}=   \begin{bmatrix}\frac{\partial N}{\partial \bar u} & \frac{\partial N}{\partial \bar \theta _{_1}} & \frac{\partial N}{\partial \bar \theta _{_2}}     \\     \frac{\partial M_{_1}}{\partial \bar u} & \frac{\partial M_{_1}}{\partial \bar \theta _{_1}} &     \frac{\partial M_{_1}}{\partial \bar \theta _{_2}}     \\     \frac{\partial M_{_2}}{\partial \bar u} & \frac{\partial M_{_2}}{\partial \bar \theta _{_1}} &     \frac{\partial M_{_2}}{\partial \bar \theta _{_2}}   \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

En los próximos apartados se obtiene la matriz de rigidez elástica para los elementos de viga de Euler-Bernoulli, de Timoshenko y de Euler-Bernoulli no lineal. Para el último elemento se tiene en cuenta el acoplamiento entre los efectos de los esfuerzos de flexión y axil lo que genera coeficientes no lineales en la matriz de rigidez elástica, por esta razón se le denomina aquí con el adjetivo ''no lineal''.
<!--
==3 Elemento de viga Euler-Bernoulli==

El movimiento deformacional de este elemento es descrito por las siguientes funciones de forma

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}u=&\frac{x}{l_{_0}}\bar u     \\     v=&x\bigl(1-\frac{x}{l_{_0}}\bigr)^{^2}\bar \theta _{_1}+\frac{x^2}{l_{_0}}     \bigl(\frac{x}{l_{_0}}-1\bigr)\bar \theta _{_2}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

Es importante destacar que los desplazamientos transversales <math display="inline">v_{_1}</math> y <math display="inline">v_{_2}</math> son nulos, esto debido a que las funciones de forma están descritas en relación al sistema de referencia local en la configuración deformada del elemento conforme se muestra en la Figura 1. Por consiguiente, la curvatura y la deformación de un punto en el interior del elemento se escriben como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\kappa =&\frac{\partial ^2 v}{\partial x^2}=     \bigl(-\frac{4}{l_{_0}}+6\frac{x}{l_{_0}^{^2}}\bigr)\bar \theta _{_1}+     \bigl(-\frac{2}{l_{_0}}+6\frac{x}{l_{_0}^{^2}}\bigr)\bar \theta _{_2}     \\     \varepsilon =&\frac{\partial u}{\partial x}-\kappa y=\frac{\bar u}{l_{_0}}+y     \biggl(\bigl(\frac{4}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)\bar \theta _{_1}+     \bigl(\frac{2}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)\bar \theta _{_2}     \biggr)   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

Al aplicar una variación virtual en el vector de desplazamientos locales <math display="inline">(\delta \bar u,\delta \bar \theta _{_1},\delta \bar \theta _{_2})</math>, se tiene una variación virtual de la deformación <math display="inline">\delta \varepsilon </math> y por el princípio de los trabajos virtuales se tiene que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}V=\int \limits _V\sigma \delta \varepsilon dV=N\delta \bar u+M_{_1}\delta \bar \theta _{_1}+M_{_2}\delta \bar \theta _{_2}     \\     N=\int \limits _A\sigma dA,~~~M_{_1}=\int \limits _A\sigma ydA,~~~M_{_2}=-\int \limits _A\sigma ydA   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

Asumiendo que <math display="inline">\sigma =E\varepsilon </math> y sustituyendo la ecuación (28b) en la ecuación (29b) se obtiene el vector de esfuerzos internos que escrito como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}_{_l}=   \begin{Bmatrix}N\\     M_{_1}\\     M_{_2}   \end{Bmatrix}=   \begin{Bmatrix}\frac{EA}{l_{_o}}\bar u \\     \frac{EI}{l_{_o}}(4\bar \theta _{_1}+2\bar \theta _{_2})\\     \frac{EI}{l_{_o}}(4\bar \theta _{_2}+2\bar \theta _{_1})   \end{Bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Ahora diferenciando cada componente del vector de esfuerzos internos en la ecuación (30) respecto a los desplazamientos locales <math display="inline">(\bar u,\bar \theta _{_1},\bar \theta _{_2})</math>, se obtiene la matriz de rigidez elástica dada por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf K}_{_l}=   \begin{bmatrix}\frac{EA}{l_{_o}} & 0 & 0 \\     0 & \frac{4EI}{l_{_o}} & \frac{2EI}{l_{_o}} \\     0 & \frac{2EI}{l_{_o}} & \frac{4EI}{l_{_o}}   \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

==4 Elemento de viga Timoshenko==

En el elemento de viga de Timoshenko son utilizadas funciones de interpolaciones lineales para los desplazamientos <math display="inline">u</math>, <math display="inline">v</math> y <math display="inline">\theta </math> en el sistema de referencia local como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}u=&\frac{x}{l_{_0}}\bar u     \\     v=&0     \\     \theta =&(1-\frac{x}{l_{_0}})\bar \theta _{_1}+\frac{x}{l_{_0}}\bar \theta _{_2}   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

La curvatura <math display="inline">\kappa </math>, la deformación por corte <math display="inline">\gamma </math> y la deformación axial <math display="inline">\varepsilon </math> son dadas por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\kappa =&\frac{\partial \theta }{\partial x}=\frac{\bar \theta _{_2}-\bar \theta _{_1}}{l_{0}}     \\     \gamma =&\frac{\partial v}{\partial x}-\theta=-(1-\frac{x}{l_{_0}})\bar \theta _{_1}-     \frac{x}{l_{_0}}\bar \theta _{_2}     \\     \varepsilon =&\frac{\partial u}{\partial x}-\kappa y=\frac{\bar u}{l_{_0}}     -\frac{\bar \theta _{_2}-\bar \theta _{_1}}{l_{_0}}y   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

Las ecuaciones constitutivas se definen como <math display="inline">\sigma =E\varepsilon </math> and <math display="inline">\tau =G\gamma </math>. Los esfuerzos internos son calculados usando el princípio de los trabajos virtuales teniendo en cuenta las deformaciones por cortante. Así, dicho princípio puede ser escrito como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V=\int \limits _V(\sigma \delta \varepsilon{+\tau}\delta \gamma )dV=N\delta \bar u+M_{_1}\delta \bar \theta _{_1}+M_{_2}\delta \bar \theta _{_2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

Además las variaciones virtuales <math display="inline">\delta \gamma </math> y <math display="inline">\delta \varepsilon </math> se calculan a partir de las ecuaciones (33b) y (33c), respectivamente, lo que permite escribir que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int \limits _V\Bigl[\sigma \Bigl(\frac{\delta \bar u}{l_{_0}}-\frac{\delta \bar \theta _{_2}-\delta \bar \theta _{_1}}{l_{_0}}y\Bigr)-     \tau \Bigl(\bigl(1-\frac{x}{l_{_0}}\bigr)\delta \bar \theta _{_1}+\frac{x}{l_{_0}}\delta \bar \theta _{_2}\Bigr)\Bigr]dV </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

Utilizando solo un punto de gauss, <math display="inline">x=\frac{l_{_0}}{2}</math>, para evitar el bloqueo por cortante en la ecuación (35) y la ecuación (34), los esfuerzos internos son obtenidos como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}N=&\int \limits _V\frac{\sigma }{l_{_0}}dV=\int \limits _A\sigma dA     \\     M_{_1}=&     \int \limits _V\Bigl(\frac{\sigma }{l_{_0}}y-\frac{\tau }{2}\Bigr)dV=\int \limits _A\sigma ydA-\frac{l_{_0}}{2}\int \limits _A\tau dA     \\     M_{_2}=&     \int \limits _V\Bigl(-\frac{\sigma }{l_{_0}}y-\frac{\tau }{2}\Bigr)dV=-\int \limits _A\sigma ydA-\frac{l_{_0}}{2}\int \limits _A\tau dA   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Teniendo en cuenta las ecuaciones constitutivas definidas anteriormente y las ecuaciones (33b) y (33c), se obtiene la siguiente expresión para el vector de esfuerzos internos

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}_{_l}=   \begin{Bmatrix}N\\     M_{_1}\\     M_{_2}   \end{Bmatrix}=   \begin{Bmatrix}\frac{EA}{l_{_o}}\bar u \\     \frac{EI}{l_{_o}}(\bar \theta _{_1}-\bar \theta _{_2})+\frac{1}{4}GAl_{_o}(\bar \theta _{_1}+\bar \theta _{_2})     \\     \frac{EI}{l_{_o}}(\bar \theta _{_2}-\bar \theta _{_1})+\frac{1}{4}GAl_{_o}(\bar \theta _{_2}+\bar \theta _{_1})   \end{Bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

Por último, diferenciando cada componente del vector de esfuerzos internos en la ecuación (37) respecto a los desplazamiento locales <math display="inline">(\bar u,\bar \theta _{_1},\bar \theta _{_2})</math> se obtiene la matriz de rigidez elástica dada por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf K}_{_l}=   \begin{bmatrix}\frac{EA}{l_{_o}} & 0 & 0 \\     0 & \frac{EI}{l_{_o}}+\frac{1}{4}GAl_{_o} & -\frac{EI}{l_{_o}}+\frac{1}{4}GAl_{_o} \\     0 & -\frac{EI}{l_{_o}}+\frac{1}{4}GAl_{_o} & \frac{EI}{l_{_o}}+\frac{1}{4}GAl_{_o}   \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

==5 Elemento de viga Euler-Bernoulli no lineal==

La denominación ''no lineal'' se refiere al acoplamiento entre los efectos de los esfuerzos axil y de flexión, esto resulta en la obtención de coefcientes no lineales en la matriz de rigidez elástica del elemento. Por lo tanto, la deformación en un punto del elemento viene dada por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varepsilon =\varepsilon _{_f}-\kappa y=\frac{1}{l_{_0}}\int \limits _0^{l_{_0}}   \biggl(\frac{\partial u}{\partial x}+\frac{1}{2}\biggl(\frac{\partial v}{\partial x}\biggr)^{2}\biggr)dx   -\kappa y </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

donde <math display="inline">\varepsilon _{_f}</math> es una medida de distribución média de deformación axial a lo largo del elemento, esto con el objetivo de evitar el bloqueo por deformación de membrana. Se nota en el intregando de la ecuación (39) que para el cálculo de dicha deformación se incluye tanto la deformación axial, <math display="inline">\frac{\partial u}{\partial x}</math>, como la rotación de la sección transversal, <math display="inline">\frac{\partial v}{\partial x}</math>. Usando las funciones de interpolación dadas por la ecuación (27) y la definición de la curvatura dada en la ecuación (28a), la ecuación anterior se escribe como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varepsilon =\frac{\bar u}{l_{_0}}+\frac{1}{15}\bar \theta _{_1}^{^2}-\frac{1}{30}\bar \theta _{_1}\bar \theta _{_2}+   \frac{1}{15}\bar \theta _{_2}^{^2}+y     \biggl(\bigl(\frac{4}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)\bar \theta _{_1}+     \bigl(\frac{2}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)\bar \theta _{_2}     \biggr) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

cuya variación virtual se expresa de la siguiente manera

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \varepsilon =\frac{\delta \bar u}{l_{_0}}+\frac{2}{15}\bar \theta _{_1}\delta \bar \theta _{_1}   -\frac{1}{30}\bar \theta _{_2}\delta \bar \theta _{_1}-\frac{1}{30}\bar \theta _{_1}\delta \bar \theta _{_2}+   \frac{2}{15}\bar \theta _{_2}\delta \bar \theta _{_2}+y     \biggl(\bigl(\frac{4}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)\delta \bar \theta _{_1}+     \bigl(\frac{2}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)\delta \bar \theta _{_2}     \biggr) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

Al aplicar el principio de los trabajos virtuales y teniendo en cuenta la ecuación (41) se obtienen los esfuerzos internos como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}N=&\int \limits _V\frac{\sigma }{l_{_0}}dV=\int \limits _A\sigma dA     \\     M_{_1}=&\Bigl(\frac{2}{15}\bar \theta _{_1}-\frac{1}{30}\bar \theta _{_2}\Bigr)\int \limits _V\sigma dV+     \int \limits _V\sigma y\bigl(\frac{4}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)dV\left.=\right.     \Bigl(\frac{2}{15}\bar \theta _{_1}-\frac{1}{30}\bar \theta _{_2}\Bigr)l_{_0}\int \limits _A\sigma dA+\int \limits _A\sigma ydA     \\     M_{_2}=&\Bigl(\frac{2}{15}\bar \theta _{_2}-\frac{1}{30}\bar \theta _{_1}\Bigr)\int \limits _V\sigma dV+     \int \limits _V\sigma y\bigl(\frac{2}{l_{_0}}-6\frac{x}{l_{_0}^{^2}}\bigr)dV\left.=\right.     \Bigl(\frac{2}{15}\bar \theta _{_2}-\frac{1}{30}\bar \theta _{_1}\Bigr)l_{_0}\int \limits _A\sigma dA-\int \limits _A\sigma ydA   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

Suponiendo que <math display="inline">\sigma =E\varepsilon </math> y utilizando la ecuación (40) se obtiene el vector de esfuerzos internos como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\bf f}_{_l}=   \begin{Bmatrix}N\\     M_{_1}\\     M_{_2}   \end{Bmatrix}=   \begin{Bmatrix}EAe \\     EAl_{_o}ee_{_1}+\frac{EI}{l_{_o}}(4\bar \theta _{_1}+2\bar \theta _{_2})\\     EAl_{_o}ee_{_2}+\frac{EI}{l_{_o}}(4\bar \theta _{_2}+2\bar \theta _{_1})   \end{Bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}e=&\frac{\bar u}{l_{_o}}+\frac{1}{15}\bar \theta _{_1}^{^2}     -\frac{1}{30}\bar \theta _{_1}\bar \theta _{_2}+\frac{1}{15}     \bar \theta _{_2}^{^2}\\     e_{_1}=&\frac{2}{15}\bar \theta _{_1}-\frac{1}{30}\bar \theta _{_2}\\     e_{_2}=&\frac{2}{15}\bar \theta _{_2}-\frac{1}{30}\bar \theta _{_1}   \end{array}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

Por último, diferenciando cada componente del vector de esfuerzos internos en la ecuación (37) respecto a los desplazamiento locales <math display="inline">(\bar u,\bar \theta _{_1},\bar \theta _{_2})</math> se obtiene la matriz de rigidez elástica dada por la siguiente expresión

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{K}_{_l}=   \begin{bmatrix}\frac{EA}{l_{_o}} & EAe_{_1} & EAe_{_2} \\     EAe_{_1} & EAl_{_o}(\frac{2}{15}e+e_{_1}^2)+\frac{4EI}{l_{_o}} &     EAl_{_o}(-\frac{1}{30}e+e_{_1}e_{_2})+\frac{2EI}{l_{_o}} \\     EAe_{_2} & EAl_{_o}(-\frac{1}{30}e+e_{_1}e_{_2})+\frac{2EI}{l_{_o}} &     EAl_{_o}(\frac{2}{15}e+e_{_2}^2)+\frac{4EI}{l_{_o}}   \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

Se puede observar que los términos <math display="inline">e</math>, <math display="inline">e_{_1}e_{_2}</math>, <math display="inline">e_{_1}^2</math> y <math display="inline">e_{_2}^2</math> son funciones cuadráticas de las variables <math display="inline">\bar \theta _{_1}</math> y <math display="inline">\bar \theta _2</math>, por lo tanto los coeficientes de la matriz en (45) que poseen estos términos son no lineales respecto a <math display="inline">\bar \theta _{_1}</math> y <math display="inline">\bar \theta _{_2}</math>.

-->
==6 Ejemplos numéricos==

Es importante resaltar que entre los tres elementos de viga presentados en este trabajo, el elemento de viga de Timoshenko es el menos refinado, pues utiliza funciones de interpolaciones lineales, mientras que el elemento de viga de Euler-Bernoulli no lineal es el más refinado pues además de utilizar funciones de interporlación de grado 3 tiene en cuenta el acoplamiento entre los efectos de los esfuerzos axil y de flexión. En consecuencia, la respuesta no lineal de un sistema estrutural obtenida utilizando los tres tipos de elementos, será la misma si las magnitudes de las deformaciones fueran inifinitesimales, mientras que si las deformaciones son pequeñas o medianas la respuesta será distinta. Para realizar el análisis no lineal geométrico de los ejemplos presentados en este apartado, se utilizó un programa escrito en lenguaje Fortran90 denominado 2Dbeam<math display="inline">\_</math>nl.f90 escrito por los autores de este artículo. Se utilizará a continuación la denominaco EB para el elemento de viga de Euler-Bernoulli, EBNL para el elemento de viga Euler-Bernoulli no lineal y TI para el elemento de viga de Timoshenko. En todos los ejemplos presentados a continuación fue adoptada una tolerancia para convergencia de <math display="inline">10^{^{-5}}</math>. El objetivo es mostrar el desempeño de estos elementos de viga para reproducir el comportamiento fuertemente no lineal de algunos sistemas estructurales. Cabe destacar, como será visto en los siguientes ejemplos, que los elementos de viga aquí descritos poseen habilidad de tratar grandes rotaciones de cuerpo rígido, y además que el elemento de viga de Timoshenko no presenta bloqueo por cortante.

===6.1 Columna de Euler===

En este ejemplo se considera el pandeo de la columna de Euler con una extremidad empotrada y la otra libre conforme la Figura 2a. En dicha figura se detallan las propiedades geométricas y mecánicas de la columna que fue discretizada con 10 elementos de viga EB, EBNL y TI, respectivamente. <div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_709676212-euler_column.png|326px|Columna de Euler. (a) Propiedades geométricas y mecánicas. (b) Trayectoria secundaria de equilibrio.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 2:''' Columna de Euler. (a) Propiedades geométricas y mecánicas. (b) Trayectoria secundaria de equilibrio.
|}

Para obtener la trayectoria secundaria de equilibrio fue aplicada en la extremidad livre de la columna una carga de <math display="inline">0.001P</math> en el sentido transversal como se muestra en la Figura 2a. En la Figura 2b se muestra dicha trayectoria obtenida por la aplicación de la integral elíptica completa de primera especie [8] que fue comparada con los caminos de equilibrio obtenidos con los tres tipos de elementos de viga. Se puede observar la buena concordancia entre los resultados analítico y numéricos, y también se observa una coincidencia entre los resultados numéricos obtenidos con los diferentes elementos de viga, lo que implica que las deformaciones son infinitesimales en el comportamiento postcrítico de la columna, aunque la misma sufre grandes traslaciones y rotaciones de cuerpo rígido. Para obtener la respuesta no lineal del pandeo de la columna de Euler se utilizó el método de control variable de los desplazamientos presentados en la referencia [9]. El valor de la longitud de arco fue de 0.01 para 500 pasos de carga. El número medio de iteraciones fue de 2.0. Estos datos fueran constantes para los tres modelos de elementos de viga utilizados en el análisis numérico.

===6.2 Pórtico de Lee===

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_709676212-lee_frame_second_example.png|386px|Pórtico de Lee. (a) Propiedades geométricas y mecánicas. (b) Trayectorias de equilibrio.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 3:''' Pórtico de Lee. (a) Propiedades geométricas y mecánicas. (b) Trayectorias de equilibrio.
|}

Este ejemplo consiste en la unión de una viga y una columna formando entre si un ángulo recto. Uno de los extremos posee un apoyo con dos restricciones mientras que el otro extremo presenta un apoyo con una restricción que permite el desplazamiento en la dirección del eje <math display="inline">y</math>. Las condiciones de carga, así como las propiedades geométricas y mecánicas están descritas en la Figura 3a. Este pórtico fue discretizado con 20 elementos de viga EB, EBNL y TI, respectivamente. En la Figura 3b se muestran las trayectorias de equilibrio, donde se describen la evolución de los desplazamientos verticales <math display="inline">v_{_1}</math> del nodo rígido del pórtico y <math display="inline">v_{_2}</math> del apoyo de una restricción, durante el proceso de carga para los tres modelos de elementos de viga. Se puede observar que en la mayor parte del camino de equilibrio las trayectorias son coincidentes, tanto para <math display="inline">v_{_1}</math> como para <math display="inline">v_{_2}</math>, excepto a los valores próximos a los puntos limites PL2, PL4 y PL5 como se muestra en Figura 3b. Las diferencias se dan debido a que las deformaciones en estos tramos de la trayectoria son medianas. Los resultados de ambas trayectorias de equilibrio fueron comparados con los obtenidos por Fuji et al. [10] donde fue usada una malla con 10 elementos de viga 2D. En la Figura 3b se muestra que hay una buena concordancia con las trayectorias obtenidas por dichos autores con excepción de los puntos limites PL2 y PL3 donde se forman unos ''loopings'' cuando se utiliza la malla de 20 elementos. Para obtener la respuesta no lineal del pórtico de Lee se utilizó el método de longitud de arco con restricción cilíndrica. La longitud de arco inicial utilizada para los tres tipos de elemento de viga fue de 12.5. Para el elemento EB fueron 285 pasos de carga con 43 cortes automáticos debido a la divergencia del proceso iterativo. El número medio de iteraciones fue de 4.43. Con el elemento EBNL fueron 300 pasos de carga, 55 cortes automáticos y un número medio de iteraciones de 4.41. Por último, para el elemento TI fueron 350 pasos de carga con 78 cortes automáticos y un número medio de iteraciones de 4.65. Se puede notar que el comportamiento del Pórtico Lee es fuertemente no lineal con la presencia de algunos ''limit points'', ''turning points'' y ''loopings'' como se muestra en la Figura 3b.

===6.3 Arco circular biarticulado rebajado===

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_709676212-hinged_low_rise_arch.png|357px|Arco circular biarticulado rebajado. (a) Propiedades geométricas y mecánicas. (b) Desplazamiento vertical del vértice.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 4:''' Arco circular biarticulado rebajado. (a) Propiedades geométricas y mecánicas. (b) Desplazamiento vertical del vértice.
|}

En este ejemplo se considera un arco rebajado de directriz circular, ambos extremos apoyados y sometido a una carga <math display="inline">q</math> uniformemente distribuida en la mitad de la luz del arco conforme se muestra en la Figura 4a. En dicha figura se detallan las propiedades geométricas y mecánicas del arco que fue discretizado con 20 elementos de viga EB, EBNL y TI, respectivamente. En la Figura 4b se muestran las curvas de la carga normalizada <math display="inline">\frac{qr^3}{EI}</math> versus el desplazamiento vertical del vértice normalizado <math display="inline">\frac{v}{r}</math>, para los tres modelos de elementos de viga. Se puede observar que no hay diferencias significativas entre las tres trayectorias de equilibrio, excepto en los valores próximos a los puntos limite PL2, PL3 y PL4, respectivamente. Los valores normalizados de la carga en los puntos limites PL1, PL2 y PL3 fueron comparados con los valores obtenidos por Xu y Mirmiran [11]. Estos autores utilizaran la formulación corrotacional y un elemento de viga que incluye los términos no lineales de la deformación axial de Green. Además, utilizaron una malla con 20 elementos. Los valores de estas cargas son detalladas en la Tabla 1. Para el elemento EB la mayor diferencia es de 3.1%, observada en el punto limite PL2. Para el elemento EBNL la mayor diferencia fue de 0.53% en el punto limite PL3. Se puede observar que para este elemento fueran encontradas las menores diferencias, esto debido a la inclusión del término <math display="inline">\varepsilon _{_f}</math> definido en la ecuación (39), en la deformación infinitesimal del elemento. Para el elemento TI la mayor diferencia fue de 5.92% en el punto limite PL2. Por otra parte, las mayores diferencias fueron detectadas para este elemento pues utiliza funciones de interpolaciones lineales para evaluar la deformación axial.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size: 75%;"> '''Tabla 1'''. Valores de la carga <math>\frac{qr^3}{EI}</math> en los puntos limites</div>

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;"
|- style="border-top: 2px solid;"
|                               
| PL1 
| PL2 
| PL3        
|- style="border-top: 2px solid;"
|       Xu y Mirmiran         
| 13.77 
| -20.09 
| 33.99 
|- style="border-top: 2px solid;"
|       EB              
| 13.92 
| -20.71 
| 34.86 
|- style="border-top: 2px solid;"
|       diferencia (%) 
| 1.09  
| 3.09   
| 2.56  
|- style="border-top: 2px solid;"
|       EBNL            
| 13.83 
| -20.25 
| 34.17 
|- style="border-top: 2px solid;"
|       diferencia (%) 
| 0.44  
| 0.8    
| 0.53  
|- style="border-top: 2px solid;"
|       TI              
| 14.02 
| -21.28 
| 35.62 
|- style="border-top: 2px solid;border-bottom: 2px solid;"
|       diferencia (%) 
| 1.82  
| 5.92   
| 4.8   

|}

Para obtener las trayectorias de equilibrio fue utilizado el método de control variable de los desplazamientos. El valor inicial de la longitud de arco fue de 0.1 para los tres tipos de elementos de viga. Para el elemento EB fue utilizado 177 pasos de carga con 12 cortes automáticos debido a la divergencia del proceso iterativo. El número medio de iteraciones fue de 4.18. Con el elemento EBNL fueron 173 pasos de carga, 10 cortes automáticos y un número medio de iteraciones de 4.20. Por último, para el elemento TI fueron 173 pasos de carga con 9 cortes automáticos y un número medio de iteraciones de 4.22. Se puede notar que las trayectorias de equilibrio del arco mostrado en la Figura 4b son fuertemente no lineales con la presencia de algunos ''limit points'', ''turning points'' y un ''looping''.

===6.4 Arco semicircular biarticulado===

En este ejemplo se considera un arco semicircular con ambos extremos apoyados, sometido a la acción de una carga <math display="inline">P</math> excéntrica respecto al vértice del arco, como se muestra en la Figura 5a. En dicha figura se detallan las propiedades geométricas y mecánicas del arco que fue discretizado con 50 elementos de viga EB, EBNL y TI. En la Figura 5b se muestran las curvas carga <math display="inline">P</math> versus el desplazamiento <math display="inline">v</math> del vértice, obtenidas para los tres modelos de elementos de viga. Se puede observar que hay una buena concordancia entre las tres trayectorias de equilibrio, excepto en los valores próximos a los puntos limite PL6, PL7, PL8  y PL9, respectivamente. Sin embargo, las diferencias se acentuan en los puntos limite PL8 y PL9, debido probablemente, a que en estos tramos de equilibrio la elástica del arco presenta deformaciones medianas. Los valores de la carga <math display="inline">P</math> en los nueve puntos limites fueron comparados con los valores presentados por Yang y Kuo [12]. Estos autores utilizaran la formulación Lagrageana Actualizada, un elemento de viga que incluye los términos no lineales de la deformación axial de Green y una malla con 26 elementos. Los valores de estas cargas se detallan en la Tabla 2. Para el elemento EB la mayor diferencia es de 3.82% observada en el punto limite LP9. Para el elemento EBNL la mayor diferencia fue de 1.76% en el punto limite PL8. Para el elemento TI la mayor diferencia ha sido de 6.88% en el punto limite PL9. Para los puntos PL1, PL2, PL3 y PL4 las diferencias no son significativas, sin embargo para los puntos PL6, PL7, PL8 y PL9 se han detectado las mayores diferencias.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_709676212-hinged_semi_circular_arch.png|427px|Arco semicircular biarticulado. a) Propiedades geométricas y mecánicas. b) Desplazamiento vertical del vértice.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 5:''' Arco semicircular biarticulado. (a) Propiedades geométricas y mecánicas. (b) Desplazamiento vertical del vértice.
|}

Para obtener las trayectorias de equilibrio fue utilizado el método de control variable de los desplazamientos. El valor inicial de la longitud de arco fue de 2.2 para los tres tipos de elementos de viga. Para el elemento EB se ha utilizado 859 pasos de carga con 91 cortes automáticos debido a la divergencia del proceso iterativo. El número medio de iteraciones fue de 4.23. Con el elemento EBNL fueron 820 pasos de carga, 93 cortes automáticos y un número medio de iteraciones de 4.25. Por último, para el elemento TI fueron 859 pasos de carga con 89 cortes automáticos y un número medio de iteraciones de 4.22. Se puede notar que las trayectorias de equilibrio del arco mostrado en la Figura 5b son fuertemente no lineales con la presencia de varios ''limit points'', ''turning points'' y ''loopings''.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;"
|-style="text-align: center; font-size: 75%;"
| colspan='10' | '''Tabla 2'''. Valores de la carga <math>P</math> en los puntos limites - (lb)

|- style="border-top: 2px solid;font-size: 85%;"
|                       
| PL1 
| PL2 
| PL3 
| PL4 
| PL5 
| PL6 
| PL7 
| PL8 
| PL9 
|- style="border-top: 2px solid;font-size: 85%;"
|      Yang y Kuo     
| 5.813 
| -8.498 
| 16.149 
| -22.162 
| 38.566 
| -49.896 
| 64.875                          
| -82.420
| 104.611
|- style="border-top: 2px solid;font-size: 85%;"
|      EB             
| 5.811 
| -8.495 
| 16.204 
| -22.086 
| 38.932 
| -50.206 
| 66.786                          
| -83.138
| 108.61 
|- style="border-top: 2px solid;font-size: 85%;"
|      diferencia (%)
| 0.03  
| 0.04   
| 0.34   
| 0.34    
| 0.95   
| 0.62    
| 2.95                          
| 0.87  
| 3.82   
|- style="border-top: 2px solid;font-size: 85%;"
|      EBNL           
| 5.802 
| -8.464 
| 16.108 
| -21.912 
| 38.453 
| -49.394 
| 65.274                          
| -80.967
| 104.99
|- style="border-top: 2px solid;font-size: 85%;"
|      diferencia (%)
| 0.19  
| 0.40   
| 0.25   
| 1.13    
| 0.29   
| 1.00    
| 0.62                          
| 1.76  
| 0.36   
|- style="border-top: 2px solid;font-size: 85%;"
|      TI             
| 5.816 
| -8.518 
| 16.278 
| -22.24  
| 39.328 
| -50.909 
| 68.081                          
| -85.055
| 111.81
|- style="border-top: 2px solid;border-bottom: 2px solid;font-size: 85%;"
|      diferencia (%)
| 0.05   
| 0.24  
| 0.80   
| 0.35    
| 1.98   
| 2.03    
| 4.94                          
| 3.20   
| 6.88  

|}

==7 Conclusiones==

Los elementos de viga 2D de Euler-Bernoulli, de Euler-Bernoulli no lineal y de Timoshenko presentaron un excelente desempeño en el análisis no lineal de vigas, pórticos y arcos. Estos elementos obtuvieron trayectorias de equilibrio casi idénticas para las diferentes tipologías estructurales presentadas en este trabajo, cuando estas sufren deformaciones infinitesimales. Cuando la estructura analisada ha presentado en algunos estados de carga deformaciones medianas, los caminos de equilibrio difieren un poco. Asi que de manera general los resultados obtenidos con los tres modelos de vigas han tenido una buena coincidencia comparados con resultados de otros autores. En este trabajo se mostró, a través de los ejemplos numéricos, la habilidad de los tres modelos viga de tratar con grandes rotaciones de cuerpo rígido cuando se utiliza la formulación corrotacional. A través de operaciones algebraicas bastante simples se describió de manera explícita el movimiento de sólido rígido del elemento de viga y la obtención de la matriz de rigidez tangente expresada en coordenadas globales. La formulación corrotacional aquí descrita es independiente del modelo de elemento de viga 2D que será utilizado. Cabe destacar como principal ventaja de la formulación corrotacional el desacoplamiento entre los efectos locales y globales, esto permite utilizar una biblioteca de elementos finitos oriundos del análisis lineal. La extensión de estos al análisis no lineal geométrico se da al tener en cuenta los efectos globales debido al movimiento de sólido rígido. Por otro lado, diferentes no linealidades físicas pueden ser fácilmente incorporadas a la formulación corrotacional debido a que son efectos locales.

==Referencias==
<div class="auto" style="width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

[1]&nbsp;C.A. Felippa and B. Haugen, A unified formulation of small-strain corotational finite elements: I. Theory. Computer Methods in Applied Mechanics and Engineering, 81:131-150, (2005).

[2]&nbsp;L.L. Yaw, N. Sukumar and  S.K. Kunnath, Meshfree co-rotational formulation for two-dimensional continua. International Journal for Numerical Methods in Engineering, 79(8):979–1003, (2009).

[3]&nbsp;Xu, J., C. K. Lee, C.K., Tan, K.H., A two-dimensional co-rotational Timoshenko beam element with XFEM formulation. Computational Mechanics,  49(5):667–683, (2012).

[4]&nbsp;Caselli, F. and Bisegna, P., A corotational flat triangular element for large strain analysis of thin shells with application to soft biological tissues. Computational Mechanics,  54(3):847–864, (2014).

[5]&nbsp;Li, Z.X., Xiang, Y., Izzuddin, B.A., Vu-Quoc, L., Zhuo, X., Zhang, C., A 6-node co-rotational triangular elasto-plastic shell element. Computational Mechanics, 55(5):837–859, (2015).

[6]&nbsp;Crisfield, M.A., Non-linear finite element analysis of solids and structures. Volume 1: Essentials, Wiley, (1991).

[7]&nbsp;S. Krenk, Non-linear modeling and analysis of solids and structures. Cambridge University Press, (2009).

[8]&nbsp;Timoshenko, S.P. y Gere, J.M., Theory of elastic stability. Second Edition, Dover, (2009).

[9]&nbsp;W.T. Matias, El control variable de los desplazamientos en el análisis no lineal elástico de estructuras de barras. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, 18(4):549–572, (2002).

[10]&nbsp;F. Fujii, K.K. Choong and S.X. Gong, Variable displacement control to overcome turnning points of nonlinear elastic frames. Computers & Structures, 44(1/2):133-136, (1992).

[11]&nbsp;Xu, Z. and Mirmiran, A., Looping behavior of arches using corotational finite element. Computers & Structures, 62(6):1059-1071, (1997).

[12]&nbsp;Yang, Y.B. and Kuo, S.R., Theory & analysis of nonlinear framed strucutures. Prentice Hall, (1994).
</div>