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

${\begin{array}{rl}{\bf {r}}=&{\bigl (}-{cos}\beta ,-{sen}\beta ,0,{cos}\beta ,{sen}\beta ,0{\bigr )}^{^{T}}\\{\bf {z}}=&{\bigl (}{sen}\beta ,-{cos}\beta ,0,-{sen}\beta ,{cos}\beta ,0{\bigr )}^{^{T}}\end{array}}$

(18 error)

Movimiento del elemento de viga 2D.

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 ${\textstyle (x_{_{1}},y_{_{1}})}$ y ${\textstyle (x_{_{2}},y_{_{2}})}$, respectivamente. En la configuración deformada dichas coordenadas se definem como ${\textstyle {(x'}_{_{1}}{,y'}_{_{1}})}$ y ${\textstyle {(x'}_{_{2}}{,y'}_{_{2}}).}$ A cada elemento se le asigna el sistema de referencia local ${\textstyle (x_{_{l}},y_{_{l}})}$ 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

${\bf {u}}_{_{g}}=(u_{_{1}},v_{_{1}},\theta _{_{1}},u_{_{2}},v_{_{2}},\theta _{_{2}})^{^{T}}$

(1)

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

${\bf {u}}_{_{l}}=({\bar {u}},{\bar {\theta }}_{_{1}},{\bar {\theta }}_{_{2}})^{^{T}}$

(2)

donde ${\textstyle {\bar {u}}}$ es el desplazamiento relativo entre los nodos 1 y 2 en la dirección del eje local ${\textstyle x_{_{l}}}$. ${\textstyle {\bar {\theta }}_{_{1}}}$ y ${\textstyle {\bar {\theta }}_{_{2}}}$ son las rotaciones deformacionales de los nodos 1 y 2, respectivamente. Las componentes de ${\textstyle {\bf {u}}_{_{l}}}$ son definidas de la siguiente manera

${\bar {u}}=l-l_{_{0}}$ ${\bar {\theta }}_{_{1}}=\theta _{_{1}}-\alpha$ ${\bar {\theta }}_{_{2}}=\theta _{_{2}}-\alpha$

(3)

donde ${\textstyle l_{_{0}}}$ y ${\textstyle l}$ denotan la longitud inicial y actual del elemento. Dichas longitudes son dadas por

$l_{_{0}}={\sqrt {x_{_{21}}^{2}+y_{_{21}}^{2}}}$ $l={\sqrt {(x_{_{21}}+u_{_{21}})^{2}+(y_{_{21}}+v_{_{21}})^{2}}}$

(4)

con

$x_{_{21}}=x_{_{2}}-x_{_{1}}$ $y_{_{21}}=y_{_{2}}-y_{_{1}}$ $u_{_{21}}=u_{_{2}}-u_{_{1}}$ $v_{_{21}}=v_{_{2}}-v_{_{1}}$

(5)

y ${\textstyle \alpha =\beta -\beta _{_{0}}}$ denota la rotación de cuerpo rígido que puede ser calculada como

${sen}\alpha ={cos}\beta _{_{0}}{sen}\beta {-}{sen}\beta _{_{0}}{cos}\beta$ ${cos}\alpha ={cos}\beta _{_{0}}{cos}\beta {+}{sen}\beta _{_{0}}{sen}\beta$

(6)

donde ${\textstyle \beta _{_{0}}}$ es la inclinación del elemento en la configuración de referencia y ${\textstyle \beta }$ 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

${cos}\beta _{_{0}}={\frac {x_{_{21}}}{l_{_{0}}}}$ ${sen}\beta _{_{0}}={\frac {y_{_{21}}}{l_{_{0}}}}$ ${cos}\beta ={\frac {x_{_{21}}+u_{_{21}}}{l}}={\frac {{x'}_{_{21}}}{l}}$ ${sen}\beta ={\frac {y_{_{21}}+v_{_{21}}}{l}}={\frac {{y'}_{_{21}}}{l}}$

(7)

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

${\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}}$

(8)

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

$\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 )}$

(9)

Esta expresión es singular si ${\textstyle y_{_{21}}^{\ast }=0}$, por lo que será asignado el valor de ${\textstyle \alpha =0}$ para ${\textstyle x_{_{21}}^{\ast }=l_{_{0}}l}$ o ${\textstyle \alpha =\pi }$ para ${\textstyle x_{_{21}}^{\ast }=-l_{_{0}}l}$. 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:

${\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}}$

(10)

La función módulo asigna un valor arbitrario en el intervalo ${\textstyle [0,2\pi [}$, 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 ${\textstyle \pm \pi ,\pm 2\pi ,...}$

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

$\delta {\bar {u}}=\delta l={\bigl (}-{cos}\beta ,-{sen}\beta ,0,{cos}\beta ,{sen}\beta ,0{\bigr )}^{^{T}}\delta {\bf {u}}_{_{g}}$ $\delta {\bar {\theta }}_{_{1}}=\delta \theta _{_{1}}-\delta \alpha =\delta \theta _{_{1}}-\delta \beta$ $\delta {\bar {\theta }}_{_{2}}=\delta \theta _{_{2}}-\delta \alpha =\delta \theta _{_{2}}-\delta \beta$

(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 ${\textstyle \beta }$ es calculada diferenciando la ecuación (3.7d) respecto a los desplazamientos globales y se escribe como

$\delta \beta ={\frac {1}{l}}{\bigl (}{sen}\beta ,-{cos}\beta ,0,-{sen}\beta ,{cos}\beta ,0{\bigr )}^{^{T}}\delta {\bf {u}}_{_{g}}$

(12)

Reagrupando en forma matricial las relaciones entre las variaciones virtuales de las variables locales y globales, se escribe la matriz de transformación ${\textstyle {\bf {B}}}$ como

$\delta {\bf {u}}_{_{l}}={\bf {B}}\delta {\bf {u}}_{_{g}}$

con

${\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}}$

(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

$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}}$

(14)

La igualdad descrita en la ecuación (14) es válida para cualquier valor arbitrario de ${\textstyle \delta {\bf {u}}_{_{g}}}$, de modo que

${\bf {f}}_{_{g}}={\bf {B}}^{^{T}}{\bf {f}}_{_{l}}$

(15)

donde el vector de esfuerzos internos ${\textstyle {\bf {f}}_{_{l}}=(N,M_{_{1}},M_{_{2}})^{^{T}}}$ 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

$\delta {\bf {f}}_{_{g}}={\bf {K}}_{_{g}}\delta {\bf {u}}_{_{g}}$

(16)

donde ${\textstyle {\bf {K}}_{_{g}}}$ es la matriz de rigidez tangente global que se obtiene diferenciando la ecuación (15) como

$\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}}$

(17)

donde ${\textstyle {\bf {b}}_{_{1}}}$, ${\textstyle {\bf {b}}_{_{2}}}$ y ${\textstyle {\bf {b}}_{_{3}}}$ son la primera, segunda y tercera columna de la matriz ${\textstyle {\bf {B}}^{^{T}}}$. Para auxiliar en el desarrollo algebraico son definidos los siguientes vectores

$r=(-{cos}\beta ,-{sen}\beta ,0,{cos}\beta ,{sen}\beta ,0)^{T}$ $z=({sen}\beta ,-{cos}\beta ,0,-{sen}\beta ,{cos}\beta ,0)^{T}$

(18)

cuyas variaciones virtuales son

$\delta {\bf {r}}={\bf {z}}\delta \beta$ $\delta {\bf {z}}=-{\bf {r}}\delta \beta$

(19)

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

$\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}}$

(20)

Con esta notación los vectores ${\textstyle {\bf {b}}_{_{1}}}$, ${\textstyle {\bf {b}}_{_{2}}}$ y ${\textstyle {\bf {b}}_{_{3}}}$ se escriben como

${\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}}$

(21)

cuyas variaciones virtuales son dadas por

$\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}}$

(22)

El primer término de la ecuación (17) es calculado al introducir la matriz de rigidez elástica local ${\textstyle {\bf {K}}_{_{l}}}$, 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

$\delta {\bf {f}}_{_{l}}={\bf {K}}_{_{l}}\delta {\bf {u}}_{_{l}}={\bf {K}}_{_{l}}{\bf {B}}\delta {\bf {u}}_{_{g}}$

(23)

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

${\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}})$

(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 ${\textstyle Q={\frac {M_{_{1}}+M_{_{2}}}{l}}}$, la matriz de rigidez geométrica global puede ser escrita mediante la siguiente expresión

${\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}})$

(25)

Por otra parte, la matriz de rigidez elastica ${\textstyle {\bf {K}}_{_{l}}}$ es obtenida por medio de la deferenciación de los esfuerzos internos ${\textstyle {\bf {f}}_{_{l}}}$ respecto a los desplazamientos locales ${\textstyle {\bf {u}}_{_{l}}}$, tal que

${\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}}$

(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

$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}}$

(27)

Es importante destacar que los desplazamientos transversales ${\textstyle v_{_{1}}}$ y ${\textstyle v_{_{2}}}$ 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

$\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 )}$

(28)

Al aplicar una variación virtual en el vector de desplazamientos locales ${\textstyle (\delta {\bar {u}},\delta {\bar {\theta }}_{_{1}},\delta {\bar {\theta }}_{_{2}})}$, se tiene una variación virtual de la deformación ${\textstyle \delta \varepsilon }$ y por el princípio de los trabajos virtuales se tiene que

${\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}}$

(29)

Asumiendo que ${\textstyle \sigma =E\varepsilon }$ y sustituyendo la ecuación (28b) en la ecuación (29b) se obtiene el vector de esfuerzos internos que escrito como

${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}}$

(30)

Ahora diferenciando cada componente del vector de esfuerzos internos en la ecuación (30) respecto a los desplazamientos locales ${\textstyle ({\bar {u}},{\bar {\theta }}_{_{1}},{\bar {\theta }}_{_{2}})}$, se obtiene la matriz de rigidez elástica dada por

${\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}}$

(31)

4 Elemento de viga Timoshenko

En el elemento de viga de Timoshenko son utilizadas funciones de interpolaciones lineales para los desplazamientos ${\textstyle u}$, ${\textstyle v}$ y ${\textstyle \theta }$ en el sistema de referencia local como

$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}}$

(32)

La curvatura ${\textstyle \kappa }$, la deformación por corte ${\textstyle \gamma }$ y la deformación axial ${\textstyle \varepsilon }$ son dadas por

$\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$

(33)

Las ecuaciones constitutivas se definen como ${\textstyle \sigma =E\varepsilon }$ and ${\textstyle \tau =G\gamma }$. 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

$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}}$

(34)

Además las variaciones virtuales ${\textstyle \delta \gamma }$ y ${\textstyle \delta \varepsilon }$ se calculan a partir de las ecuaciones (33b) y (33c), respectivamente, lo que permite escribir que

$\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$

(35)

Utilizando solo un punto de gauss, ${\textstyle x={\frac {l_{_{0}}}{2}}}$, para evitar el bloqueo por cortante en la ecuación (35) y la ecuación (34), los esfuerzos internos son obtenidos como

$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$

(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

${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}}$

(37)

Por último, diferenciando cada componente del vector de esfuerzos internos en la ecuación (37) respecto a los desplazamiento locales ${\textstyle ({\bar {u}},{\bar {\theta }}_{_{1}},{\bar {\theta }}_{_{2}})}$ se obtiene la matriz de rigidez elástica dada por

${\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}}$

(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

$\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$

(39)

donde ${\textstyle \varepsilon _{_{f}}}$ 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, ${\textstyle {\frac {\partial u}{\partial x}}}$, como la rotación de la sección transversal, ${\textstyle {\frac {\partial v}{\partial x}}}$. 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

$\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 )}$

(40)

cuya variación virtual se expresa de la siguiente manera

$\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 )}$

(41)

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

$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$

(42)

Suponiendo que ${\textstyle \sigma =E\varepsilon }$ y utilizando la ecuación (40) se obtiene el vector de esfuerzos internos como

${\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}}$

(43)

con

$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}}$

(44)

Por último, diferenciando cada componente del vector de esfuerzos internos en la ecuación (37) respecto a los desplazamiento locales ${\textstyle ({\bar {u}},{\bar {\theta }}_{_{1}},{\bar {\theta }}_{_{2}})}$ se obtiene la matriz de rigidez elástica dada por la siguiente expresión

${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}}$

(45)

Se puede observar que los términos ${\textstyle e}$, ${\textstyle e_{_{1}}e_{_{2}}}$, ${\textstyle e_{_{1}}^{2}}$ y ${\textstyle e_{_{2}}^{2}}$ son funciones cuadráticas de las variables ${\textstyle {\bar {\theta }}_{_{1}}}$ y ${\textstyle {\bar {\theta }}_{2}}$, por lo tanto los coeficientes de la matriz en (45) que poseen estos términos son no lineales respecto a ${\textstyle {\bar {\theta }}_{_{1}}}$ y ${\textstyle {\bar {\theta }}_{_{2}}}$.

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${\textstyle \_}$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 ${\textstyle 10^{^{-5}}}$. 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.

Columna de Euler. (a) Propiedades geométricas y mecánicas. (b) Trayectoria secundaria de equilibrio.

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 ${\textstyle 0.001P}$ 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

Pórtico de Lee. (a) Propiedades geométricas y mecánicas. (b) Trayectorias de equilibrio.

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 ${\textstyle y}$. 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 ${\textstyle v_{_{1}}}$ del nodo rígido del pórtico y ${\textstyle v_{_{2}}}$ 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 ${\textstyle v_{_{1}}}$ como para ${\textstyle v_{_{2}}}$, 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

Arco circular biarticulado rebajado. (a) Propiedades geométricas y mecánicas. (b) Desplazamiento vertical del vértice.

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 ${\textstyle q}$ 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 ${\textstyle {\frac {qr^{3}}{EI}}}$ versus el desplazamiento vertical del vértice normalizado ${\textstyle {\frac {v}{r}}}$, 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 ${\textstyle \varepsilon _{_{f}}}$ 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.

Tabla 1. Valores de la carga ${\frac {qr^{3}}{EI}}$ en los puntos limites

	PL1	PL2	PL3
Xu y Mirmiran	13.77	-20.09	33.99
EB	13.92	-20.71	34.86
diferencia (%)	1.09	3.09	2.56
EBNL	13.83	-20.25	34.17
diferencia (%)	0.44	0.8	0.53
TI	14.02	-21.28	35.62
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 ${\textstyle P}$ 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 ${\textstyle P}$ versus el desplazamiento ${\textstyle v}$ 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 ${\textstyle P}$ 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.

Arco semicircular biarticulado. a) Propiedades geométricas y mecánicas. b) Desplazamiento vertical del vértice.

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.

Tabla 2. Valores de la carga $P$ en los puntos limites - (lb)
	PL1	PL2	PL3	PL4	PL5	PL6	PL7	PL8	PL9
Yang y Kuo	5.813	-8.498	16.149	-22.162	38.566	-49.896	64.875	-82.420	104.611
EB	5.811	-8.495	16.204	-22.086	38.932	-50.206	66.786	-83.138	108.61
diferencia (%)	0.03	0.04	0.34	0.34	0.95	0.62	2.95	0.87	3.82
EBNL	5.802	-8.464	16.108	-21.912	38.453	-49.394	65.274	-80.967	104.99
diferencia (%)	0.19	0.40	0.25	1.13	0.29	1.00	0.62	1.76	0.36
TI	5.816	-8.518	16.278	-22.24	39.328	-50.909	68.081	-85.055	111.81
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

[1] 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] 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] 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] 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] 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] Crisfield, M.A., Non-linear finite element analysis of solids and structures. Volume 1: Essentials, Wiley, (1991).

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

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

[9] 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] 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] Xu, Z. and Mirmiran, A., Looping behavior of arches using corotational finite element. Computers & Structures, 62(6):1059-1071, (1997).

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