Nonlinear analysis of plane frames using a co-rotating Timoshenko beam element

Resumen

En este trabajo se describe la formulación corrotacional de un elemento de viga unificado que integra las teorías de vigas de Euler-Bernoulli y de Timoshenko y que no presenta bloqueo por deformación a cortante. La cinemática corrotacional se basa en la separación del movimiento de un sólido en una parte deformacional, y la otra, en movimiento de cuerpo rígido. La parte deformacional del movimiento es representada por tres modos de deformaciones naturales que son generados por el esfuerzo axil, la flexíon pura y la flexíon simple, respectivamente. Los esfuerzos generados por los modos naturales de deformación son auto-equilibrados lo que posibilita la obtención de la matriz de rigidez tangente corrotacional consistente. Se describe de forma detallada la obtención de las matrices de rigideces elástica, geométrica y corrotacional. A través de algunos ejemplos numéricos se muestra la habilidad del elemento de tratar con grandes rotaciones de cuerpo rígido.

Abstract

The present work describes a co-rotating shear flexible beam element without shear locking and integrating Euler-Bernoullis and Timoshenkos beam theories. The co-rotational kinematics is based on the separation of the motion in deformational and rigid body components. The deformation of the beam element is composed by three natural modes of deformation: the extension mode, the symmetric bending mode, and the anti-symmetric bending mode. The respective generalized stresses from these natural modes are self-balanced allowing the achievement of a consistent tangent stiffness matrix. In this paper, it is detailed and deduced all the algebraic steps for the deduction of the elastic stiffness matrix, the geometric stiffness matrix, and the co-rotation stiffness matrix. Some examples are presented and the numerical results demonstrate that the beam element here presented is able to handle large rotations.

Palabras clave

Elemento de viga Bernoulli/Timoshenko ; Formulación corrotacional ; Modos de deformación naturales

Keywords

Bernoulli/Timoshenko beam element ; Corotational kinematic ; Deformations natural modes

1. Introdución

En las últimas dos décadas se ha observado un creciente interés de la comunidad de mecánica computacional en la aplicación de la formulación corrotacional en el análisis no lineal de estructuras menos comunes utilizadas en sectores de tecnología de punta como en la industria aeronáutica, aeroespacial y petrolera. La idea principal de la formulación corrotacional es des componer el movimiento de un sólido en la suma de un movimiento de cuerpo rígido y un desplazamiento 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. Así, en esta formulación se asume que tanto las traslaciones cuanto 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 la no linealidad geométrica, siendo esta la principal motivación del uso de la formulación corrotacional. En los años 70 y 80 esta formulación fue utilizada en el análisis no lineal de una clase de estructuras utilizadas en la ingeniería aeroespacial (aviones, satélites, coehetes, antenas), en la ingeniería naval (navios, plataformas offshore), así como en la ingeniería mecánica (robótica, máquinas de excavaciones), teniendo poco impacto en la ingeniería civil para esta época. A mediados de los años 90 la formulación fue extendida al análisis de modelos de fallos locales y aplicaciones abarcando la dinámica no lineal. En el trabajo de Skallerud et al.[1] es utilizada la formulación corrotacional para analizar tuberías de gas submarinas expuestas a grandes cambios geométricos bien como a fallos locales. Yaw et al.[2] muestran aplicaciones de la formulación corrotacional y los métodos sin malla en el análisis de sólidos 2D incluyendo la no linealidad geométrica y física. Felippa y Haugen [3] describen una formulación corrotacional unificada para análisis estático no lineal de sólidos discretizados por medio de elementos de vigas, placas y láminas. En el trabajo de Mostafa et al.[4] es extendida la formulación para análisis estáticos utilizando elementos finitos de sólidos 2D, sólidos 3D y elementos denominados de sólidos-lámina de bajo orden con campos de deformaciones mejorados, que poseen apenas grados de libertad de traslación. Por otro lado, Matias y Bezerra [5] propusieron una formulación unificada para la descripción corrotacional de elementos de estructuras articuladas espaciales y de viga 2D. Finalmente, pueden ser citados otros autores [6]  ;  [7] que han aplicado la formulación corrotacional para análisis dinámico no lineal de sistemas estructurales discretizados con elementos de vigas 2D y 3D. En este trabajo se presenta la formulación corrotacional de un elemento de viga 2D con las siguientes particularidades: a) se define un vector de desplazamientos de cuerpo rígido del elemento a través de un sistema de referencia local fijo en su centroide; b) este elemento puede sufrir rotación de cuerpo rígido de cualquier magnitud debido a la aplicación de la función módulo; c) serán adoptados 3 modos de deformaciones naturales [8]  ;  [9] para representar la parte deformacional del movimiento del elemento; d) estos modos de deformaciones naturales representan esfuerzos internos auto-equilibrados; e) se obtiene el vector de fuerzas nodales a partir de los esfuerzos internos auto-equilibrados; f) es utilizado el principio de los trabajos virtuales complementar [10]  ;  [11] para obtener la matriz de rigidez elástica que integra las teorías de vigas de Euler-Bernoulli y de Timoshenko; g) este elemento será aquí denominado de elemento de viga 2D unificado; h) este elemento no presenta bloqueo por deformación al corte; i) se describen detalladamente el desarrollo algebraico para el cálculo de las matrices de rigidez corrotacional y geométrica, así como la expresión completa de la matriz de rigidez tangente. Finalmente, se aplica este elemento para el análisis no lineal de pórticos planos. Son presentados ejemplos numéricos en los cuales se constata la inexistencia de bloqueo por cortante en el análisis no lineal.

2. Descripción corrotacional

Sea un sistema de coordenadas globales cuya base es ortonormal conforme se muestra en la figura 1 . Para expresar las variables cinemáticas en la configuración no deformada se utilizan las coordenadas materiales (X , Y ), mientras que, en la configuración deformada las coordenadas espaciales (x , y ). En la configuración no deformada las coordenadas nodales del elemento son dadas por (X₁ , Y₁ ) y (X₂ , Y₂ ), respectivamente. Su longitud e inclinación iniciales son dadas por:

Figura 1. Movimiento del elemento de viga 2D.

$l_{0}={\sqrt {X_{21}^{2}+Y_{21}^{2}}}$ ${\varphi }_{0}=2\quad arctan({\frac {l_{0}-X_{21}}{Y_{21}}})$

( 1)

donde X₂₁  = X₂  − X₁ y Y₂₁  = Y₂  − Y₁ . Para obtener la expresión (1b) fue utilizada la relación ${\textstyle tan({\frac {1}{2}}{\varphi }_{0})}$ . Esta expresión es singular si Y₂₁  = 0, por lo que será asignado el valor de φ₀  = 0 para X₂₁  = l₀ o φ₀  = π para X₂₁  = − l₀ . De manera similar, en la configuración deformada las coordenadas nodales del elemento son dadas por (x₁ , y₁ ) y (x₂ , y₂ ), respectivamente. Su longitud e inclinación actuales son dadas por:

$l={\sqrt {x_{21}^{2}+y_{21}^{2}}}$ $\varphi =2\quad arctan({\frac {l-x_{21}}{y_{21}}})$

( 2)

donde x₂₁  = x₂  − x₁ e y₂₁  = y₂  − y₁ . Nuevamente, para obtener la expresión (2b) fue utilizada la relación ${\textstyle tan({\frac {1}{2}}\varphi }$ ). Esta expresión es singular si y₂₁  = 0, donde φ  = 0 para x₂₁  = l o φ  = π para x₂₁  = − l . Como se muestra en la figura 1 el movimiento del elemento de viga 2D es dado por el vector de desplazamientos nodales en coordenadas globales como:

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

( 3)

Por otro lado, las coordenadas espaciales nodales pueden ser escritas en función de los desplazamientos nodales como x₁  = X₁  + u₁ y x₂  = X₂  + u₂ , respectivamente. A partir de estas expresiones se obtiene que x₂₁  = X₂₁  + u₂₁ y y₂₁  = Y₂₁  + u₂₁ , donde u₂₁  = u₂  − u₁ y ${\textstyle v_{21}=v_{2}-v_{1}}$ , respectivamente, lo que permite reescribir que: ${\textstyle l={\sqrt {{\left(X_{21}+u_{21}\right)}^{2}+{\left(Y_{21}+v_{21}\right)}^{2}}}}$ .

Para definir el movimiento de cuerpo rígido de este elemento se fija en su centroide un sistema de coordenadas locales (x_e , y_e ) que lo acompaña durante todo su movimiento, como se muestra en la figura 1 . Así, el movimiento de cuerpo rígido es definido por dos traslaciones y una rotación del sistema de referencia local entre las configuraciones deformada y no deformada, y se escribe como:

${\boldsymbol {u}}_{r}=\left\{{\begin{array}{c}{\frac {1}{2}}(u_{2}+u_{1})\\{\frac {1}{2}}(v_{2}+v_{1})\\\varphi -{\varphi }_{0}\end{array}}\right\}$

( 4)

2.1. Modos de deformación naturales

Fueron necesarios tres grados de libertad para definir el movimiento de cuerpo rígido del elemento de viga 2D, así los tres grados de libertad restantes son necesarios para definir el movimiento deformacional del elemento. En este trabajo se adopta el concepto de los modos naturales de deformación definidos en las referencias [8]  ;  [9] . Para el caso del elemento de viga 2D son definidos tres modos de deformaciones naturales como se muestra en la figura 2 . El primer modo es obtenido por la aplicación de un esfuerzo axil de tracción en los nodos del elemento. El segundo modo es obtenido al considerar el elemento sujeto a un estado de flexión pura lo que implica en ausencia de esfuerzo cortante. Este modo define una línea elástica simétrica. En el tercero modo el elemento está sujeto a un estado de flexión simple que resulta en esfuerzo cortante constante. Este modo define una línea elástica anti-simétrica. Por lo tanto, se define el vector de desplazamientos deformacionales en relación al sistema de coordenadas locales como:

Figura 2. Modos naturales de deformación. a) estiramiento. b) elástica simétrica. c) elástica anti-simétrica.

${\boldsymbol {u}}_{d}^{T}=(u,{\theta }_{s},{\theta }_{a})$

( 5)

donde

$u=l-l_{0}$ ${\theta }_{s}={\theta }_{2}-{\theta }_{1}$ ${\theta }_{a}={\theta }_{2}+{\theta }_{1}-2(\varphi -{\varphi }_{0})$

( 6)

Para obtener las expresiones (6b) y (6c) se utilizan las rotaciones deformacionales de los nodos del elemento que de acuerdo con la figura 1 se definen como θ₁  − (φ  − φ₀ ) y θ₂  − (φ  − φ₀ ), respectivamente. Luego, conforme las figuras 2 b y 2 c, se suman las rotaciones simétricas y anti-simétricas de cada nodo del elemento. Por último, se igualan esas rotaciones a las rotaciones deformacionales de los nodos del elemento descritas en la figura 1 . Este procedimiento genera el siguiente sistema de ecuaciones:

$-{\frac {1}{2}}{\theta }_{s}+{\frac {1}{2}}{\theta }_{a}={\theta }_{1}-$$(\varphi -{\varphi }_{0})$ ${\frac {1}{2}}{\theta }_{s}+{\frac {1}{2}}{\theta }_{a}={\theta }_{2}-$$(\varphi -{\varphi }_{0})$

( 7)

La resolución de este sistema obtiene los valores de θ_s y θ_a descritos en las ecuaciones (6b) y (6c), respectivamente. Con el objeto de asumir cualquier magnitud en el valor de los ángulos, la expresión para rotación anti-simétrica debe ser calculada usando la función módulo [11] , tal que:

${\theta }_{a}={\underset {2\pi }{mod}}({\theta }_{a}+\pi )-$$\pi$

( 8)

La función módulo coloca un valor arbitrario en el intervalo [0, 2π [, 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 2D sufra rotaciones ±π , ± 2π , . . .

2.2. Vector de esfuerzos internos

Es importante resaltar que los modos de deformaciones naturales generan esfuerzos seccionales auto-equilibrados como se muestra en la figura 3 . Los modos de deformaciones naturales así como sus respectivos esfuerzos se expresan en relación al sistema de referencia local. Por tanto, el vector de esfuerzos auto-equilibrados es dado por:

Figura 3. Esfuerzos auto-equilibrados. a) Esfuerzo axil de tracción. b) Flexión pura. c) Flexión simple.

${\boldsymbol {f}}_{d}^{T}=(N,M_{s},M_{a})$

( 9)

La combinación de los esfuerzos seccionales auto-equilibrados, mostrados en la figura 3 , genera un vector de fuerzas nodales auto-equilibradas como se muestra en la figura 4 c. La relación de este vector con el vector de fuerzas nodales en coordenadas locales, mostrado en la figura 4 b, es dada por la siguiente expresión:

Figura 4. Vector de fuerzas nodales. a) En coordenadas globales. b) En coordenadas locales. c) Auto-equilibradas.

${\boldsymbol {f}}^{e}={\boldsymbol {S}}{\boldsymbol {f}}_{d}$

( 10)

donde S es la matriz de transformación dada por:

${\boldsymbol {S}}=\left[{\begin{array}{ccc}-1&0&0\\0&0&2/l\\0&-1&1\\1&0&0\\0&0&-2/l\\0&1&1\end{array}}\right]$

( 11)

Análogamente, la relación entre el vector de fuerzas nodales en coordenadas locales y el vector de fuerzas nodales en coordenadas globales, mostrado en la figura 4 a, es dada por:

${\boldsymbol {f}}={\boldsymbol {R}}\quad {\boldsymbol {f}}^{e}$

( 12)

donde R es la matriz de rotación del sistema de coordenadas locales para el sistema de coordenadas globales dada por:

${\boldsymbol {R}}=\left[{\begin{array}{cccccc}cos\varphi &-sen\varphi &0&0&0&0\\sen\varphi &cos\varphi &0&0&0&0\\0&0&1&0&0&0\\0&0&0&cos\varphi &-sen\varphi &0\\0&0&0&sen\varphi &cos\varphi &0\\0&0&0&0&0&1\end{array}}\right]$

( 13)

Finalmente, substituyendo la ecuación (10) en (12) , se obtiene la relación entre el vector de fuerzas auto-equilibradas y el vector de fuerzas internas en coordenadas globales que se escribe como:

${\boldsymbol {f}}={\boldsymbol {R}}\quad {\boldsymbol {S}}\quad {\boldsymbol {f}}_{d}$

( 14)

Nótese que los desplazamientos nodales pueden ser escritos en relación al sistema de coordenadas locales como:

${\boldsymbol {u}}^{e}={\boldsymbol {R}}^{T}\quad {\boldsymbol {u}}$

( 15)

3. Obtención de la matriz de rigidez tangente

Basados en la cinemática corrotacional es necesario deducir la expresión de la matriz de rigidez tangente en relación al sistema de referencia local fijo en el elemento. La matriz de rigidez tangente expresa la relación entre los incrementos de fuerzas y de desplazamientos nodales en coordenadas locales como:

$d{\boldsymbol {f}}^{e}={\boldsymbol {K}}^{e}d{\boldsymbol {u}}^{e}$

( 16)

Utilizando las ecuaciones (12) y (15) , se obtiene la matriz de rigidez tangente expresada en relación al sistema de coordenadas globales de la siguiente manera:

${\boldsymbol {K}}={\boldsymbol {R}}\quad {\boldsymbol {K}}^{e}{\boldsymbol {R}}^{T}$

( 17)

Por lo tanto, si se considera el elemento de viga 2D en su configuración deformada actual como se muestra en la figura 5 , y se aplica un incremento infinitesimal de desplazamientos a partir de esta configuración de equilibrio, la variación de la rotación de cuerpo rígido y de la longitud del elemento se escriben:

$d\varphi ={\frac {{dv}_{2}^{e}-{dv}_{1}^{e}}{l}},\quad dl={du}_{2}^{e}-$${du}_{1}^{e}$

( 18)

donde ${\textstyle ({du}_{1}^{e},{dv}_{1}^{e},{du}_{2}^{e},{dv}_{2}^{e})}$ son las componentes del incremento de desplazamientos nodales en relación al sistema de coordenadas locales. Por otro lado, recordando que el trabajo interno no varía en relación a los diferentes sistemas de coordenadas y teniendo en cuenta la ecuación (10) , la variación del trabajo interno puede ser escrita como:

Figura 5. Movimiento incremental del elemento de viga 2D.

$dV=d{\boldsymbol {u}}^{e^{T}}\quad {\boldsymbol {f}}^{e}=d{\boldsymbol {u}}^{e^{T}}{\boldsymbol {S}}\quad {\boldsymbol {f}}_{d}=$$d{\boldsymbol {u}}_{d}^{T}{\boldsymbol {f}}_{d}$

( 19)

donde ${\textstyle d{\boldsymbol {u}}_{d}^{T}=(du,d{\theta }_{s},d{\theta }_{s})}$ es la variación de los modos de deformaciones naturales. La última igualdad se cumple para cualquier valor arbitrario de f_d , en consecuencia se obtiene que:

$d{\boldsymbol {u}}_{d}={\boldsymbol {S}}^{T}d{\boldsymbol {u}}^{e}$

( 20)

La variación del vector de fuerzas internas dada por la ecuación (14) se escribe como:

$d{\boldsymbol {f}}={\boldsymbol {R}}{\boldsymbol {S}}d{\boldsymbol {f}}_{d}+$${\boldsymbol {R}}d{\boldsymbol {S}}{\boldsymbol {f}}_{d}+d{\boldsymbol {R}}{\boldsymbol {S}}{\boldsymbol {f}}_{d}$

( 21)

con

$d{\boldsymbol {R}}=\left[{\begin{array}{cccccc}-sin\varphi &-cos\varphi &0&0&0&0\\cos\varphi &-sin\varphi &0&0&0&0\\0&0&0&0&0&0\\0&0&0&-sin\varphi &-cos\varphi &0\\0&0&0&cos\varphi &-sin\varphi &0\\0&0&0&0&0&0\end{array}}\right]d\varphi ,\quad d{\boldsymbol {S}}=$$\left[{\begin{array}{ccc}0&0&0\\0&0&-2/l^{2}\\0&0&0\\0&0&0\\0&0&2/l^{2}\\0&0&0\end{array}}\right]dl$

( 22)

La relación entre las variaciones de los esfuerzos seccionales auto-equilibrados y de los modos de deformaciones naturales es dada por:

$d{\boldsymbol {f}}_{d}={\boldsymbol {K}}_{d}d{\boldsymbol {u}}_{d}$

( 23)

donde K_d es la matriz de rigidez elástica del elemento, cuya dimensión es 3 × 3. Teniendo en cuenta las ecuaciones (12) , (20) y (23) , la ecuación (21) puede ser reescrita como:

$d{\boldsymbol {f}}^{e}={\boldsymbol {S}}{\boldsymbol {K}}_{d}{\boldsymbol {S}}^{T}\quad d{\boldsymbol {u}}^{e}+$$(d{\boldsymbol {S}}+{\boldsymbol {R}}^{T}d{\boldsymbol {R}}\quad {\boldsymbol {S}}){\boldsymbol {f}}_{d}$

( 24)

Finalmente, teniendo en cuenta la ecuación (16) , la matriz de rigidez tangente en relación al sistema de coordenadas locales puede ser escrita como:

${\boldsymbol {K}}^{e}={\boldsymbol {S}}{\boldsymbol {K}}_{d}{\boldsymbol {S}}^{T}+$${\boldsymbol {K}}_{r}$

( 25)

donde el primer término tiene en cuenta los coeficientes de rigidez de los modos de deformaciones naturales, mientras que K_r es la matriz de rigidez corrotacional, que representa el efecto de la variación de la rotación de cuerpo rígido del sistema de referencia local y del efecto del esfuerzo cortante debido a la variación en la longitud del elemento, que de acuerdo con la ecuación (24) , se define como: K_r  = (d''S  + RTd''R  S )f_d . Nótese que de acuerdo con la ecuación (22) la variación de la matriz de transformación S depende de la variable dl y que la variación de la matriz de rotación R depende de dφ . Teniendo en cuenta las ecuaciones (18a) y (18b), después de algunos desarrollos algebraicos se llega a:

${\boldsymbol {K}}_{r}={\frac {1}{l}}\left[{\begin{array}{cccccc}0&Q&0&0&-Q&0\\Q&N&0&-Q&-N&0\\0&0&0&0&0&0\\0&-Q&0&0&Q&0\\-Q&-N&0&Q&N&0\\0&0&0&0&0&0\end{array}}\right]$

( 26)

Esta matriz forma parte de la matriz de rigidez geometrica completa del elemento de viga 2D como se demostrará en el apartado 3.3.

3.1. Matriz de rigidez elástica

El primer modo de deformación natural es la elongación del elemento de viga 2D y su coeficiente de rigidez, de acuerdo con la teoría de la resistencia de materiales, se da por la relación incremental ${\textstyle dN={\frac {EA}{l}}du}$ . Para el elemento en flexión pura se obtiene el coeficiente de rigidez aplicándose el Principio de Trabajos Virtuales Complementar (PTVC) que resulta en la siguiente relación incremental ${\textstyle {dM}_{s}={\frac {EI}{l}}d{\theta }_{s}}$ . El tercer modo de deformación natural representa el elemento sujeto a flexión simple que resulta en un esfuerzo cortante constante. De esta manera, al aplicarse el PTVC, se obtiene que:

$M_{a}{\theta }_{a}={\int }_{0}^{l}\left({\frac {M{\left(x\right)}^{2}}{EI}}+\right.$$\left.{\frac {Q^{2}}{{GA}_{0}}}\right)dx={\frac {1}{3}}{\frac {M_{a}^{2}l}{EI}}+$${\frac {Q^{2}l}{{GA}_{0}}}$

( 27)

donde A₀ es el área de la sección transversal corregida por el factor de forma que tiene en cuenta el efecto de la distribución de las tensiones de cortante en la sección transversal. Se introduce el efecto de la deformación por cortante a través de los siguientes parámetros:

${\psi }_{a}={\frac {1}{1+\Phi }},\quad \Phi ={\frac {12EI}{{GA}_{0}l^{2}}}$

( 28)

Y usando la expresión del esfuerzo cortante ${\textstyle Q={\frac {2M_{a}}{l}}}$ y la ecuación (27) , se determina la siguiente relación incremental:

${dM}_{a}=3{\psi }_{a}{\frac {EI}{l}}d{\theta }_{a}$

( 29)

Finalmente, se define la matriz de rigidez elástica como:

${\boldsymbol {K}}_{d}=\left[{\begin{array}{ccc}{\frac {EA}{l}}&0&0\\0&{\frac {EI}{l}}&0\\0&0&3{\psi }_{a}{\frac {EI}{l}}\end{array}}\right]$

( 30)

Cabe destacar que los coeficientes de rigidez a flexión de esta matriz tienen en cuenta tanto la hipótesis de Euler-Bernoulli cuanto la hipótesis de Timoshenko, de ahí la denominación de elemento de viga unificado.

3.2. Matriz de rigidez geométrica

Además de la matriz de rigidez elástica es necesario tener en cuenta la matriz de rigidez geométrica debido al efecto del esfuerzo axial cuando el elemento de viga 2D sufre un incremento de rotación de cuerpo rígido. En este trabajo, para obtener los coeficientes de la matriz de rigidez geométrica se considera la ecuación diferencial de la viga-columna sin cargas transversales y se desprecia el efecto a cortante. Es considerado el elemento de viga simplemente apoyado para que se cumplan las condiciones de contorno del segundo y tercer modo de deformación natural como se muestra en las figuras 2 b y 2 c. De esta manera, la ecuación diferencial de gobierno se escribe como:

${EIv}^{{''}{''}}-{Nv}^{''}=0$

( 31)

Después de aplicar dos veces la integración por partes, se obtiene la siguiente forma integral de esta ecuación:

${\int }_{0}^{l}(\delta v^{''}{EIv}^{''}+\delta v^{'}{Nv}^{'})dx=$${\left[\delta v\quad Q+\delta \theta M\right]}_{0}^{l}$

( 32)

La relación incremental de esta forma integral, después de algunos pasos algebraicos, se expresa como:

${\int }_{0}^{l}(\delta v^{''}EI\quad {dv}^{''}+\delta v^{'}{Ndv}^{'})dx=$${\left[\delta v\quad dQ+\delta \theta dM\right]}_{0}^{l}$

( 33)

Para calcular el coeficiente de rigidez geométrico debido al segundo modo de deformación natural se aplica la siguiente parte de la forma integral descrita arriba:

${\int }_{0}^{l}\delta v^{'}{Ndv}^{'}dx={\left[\delta {\theta }_{s}{dM}_{s}\right]}_{0}^{l}$

( 34)

Se interpola tanto el campo de desplazamiento incremental dv   cuanto el campo de desplazamiento virtual ${\textstyle \delta v}$ y sus derivadas de primer orden utilizando la función de forma para la elástica simétrica definida en la figura 6 . Después de aplicarse la regla de la cadena se obtiene que ${\textstyle {dv}^{'}={\frac {1}{2}}\xi d{\theta }_{s}}$ y ${\textstyle \delta v^{'}={\frac {1}{2}}\xi \delta {\theta }_{s}}$ . En seguida, integrando la ecuación (34) en el espacio de la coordenada ξ se llega a:

Figura 6. Funciones de forma.

${dM}_{a}={\frac {1}{12}}Nld{\theta }_{a}\quad \quad \Rightarrow \quad \quad k_{a}=$${\frac {1}{20}}Nl$

( 35)

Para calcular el coeficiente de rigidez geométrico debido al tercero modo de deformación natural, se aplica la siguiente parte de la forma integral descrita en la ecuación (33) :

${\int }_{0}^{l}\delta v^{'}{Ndv}^{'}dx={\left[\delta {\theta }_{a}{dM}_{a}\right]}_{0}^{l}$

( 36)

Siguiendo los pasos anteriores y recordando que para el caso de la elástica anti-simétrica al aplicar la regla de la cadena se obtiene que ${\textstyle {dv}^{'}=-{\frac {1}{4}}(1-3{\xi }^{2})d{\theta }_{a}}$ y ${\textstyle \delta v^{'}=-{\frac {1}{4}}(1-3{\xi }^{2})\delta {\theta }_{a}}$ . De esta manera, integrando la ecuación (36) en el espacio de la coordenada ξ se llega a:

${dM}_{a}={\frac {1}{12}}Nld{\theta }_{a}\quad \quad \Rightarrow \quad \quad k_{a}=$${\frac {1}{20}}Nl$

( 37)

Y se define la matriz de rigidez geométrica como:

${\boldsymbol {K}}_{g}=\left[{\begin{array}{ccc}0&0&0\\0&{\frac {1}{12}}Nl&0\\0&0&{\frac {1}{20}}Nl\end{array}}\right]$

( 38)

3.3. Matriz de rigidez tangente completa

Para obtener la matriz de rigidez tangente completa del elemento de viga 2D se extienden las matrices K_d y K_g para sus formas completas, esto es, para matrices de orden 6 × 6 utilizándose la matriz de transformación S . Por consiguiente, la forma completa se escribe como:

${\boldsymbol {K}}^{e}={\boldsymbol {S}}{\boldsymbol {K}}_{d}{\boldsymbol {S}}^{T}+$${\boldsymbol {S}}{\boldsymbol {K}}_{g}{\boldsymbol {S}}^{T}+{\boldsymbol {K}}_{r}=$${\boldsymbol {K}}_{d}^{e}+{\boldsymbol {K}}_{g}^{e}$

( 39)

con

${\boldsymbol {K}}_{d}^{e}=\left[{\begin{array}{cccccc}{\frac {EA}{l}}&0&0&-{\frac {EA}{l}}&0&0\\0&12{\psi }_{a}{\frac {EI}{l^{3}}}&6{\psi }_{a}{\frac {EI}{l^{2}}}&0&-12{\psi }_{a}{\frac {EI}{l^{3}}}&6{\psi }_{a}{\frac {EI}{l^{2}}}\\0&6{\psi }_{a}{\frac {EI}{l^{2}}}&(3{\psi }_{a}+1){\frac {EI}{l}}&0&-6{\psi }_{a}{\frac {EI}{l^{2}}}&(3{\psi }_{a}-1){\frac {EI}{l}}\\-{\frac {EA}{l}}&0&0&{\frac {EA}{l}}&0&0\\0&-12{\psi }_{a}{\frac {EI}{l^{3}}}&-6{\psi }_{a}{\frac {EI}{l^{2}}}&0&12{\psi }_{a}{\frac {EI}{l^{3}}}&-6{\psi }_{a}{\frac {EI}{l^{2}}}\\0&6{\psi }_{a}{\frac {EI}{l^{2}}}&(3{\psi }_{a}-1){\frac {EI}{l}}&0&-6{\psi }_{a}{\frac {EI}{l^{2}}}&(3{\psi }_{a}+1){\frac {EI}{l}}\end{array}}\right]$ ${\boldsymbol {K}}_{g}^{e}={\boldsymbol {S}}{\boldsymbol {K}}_{g}{\boldsymbol {S}}^{T}+$${\boldsymbol {K}}_{r}={\frac {1}{l}}\left[{\begin{array}{cccccc}0&Q&0&0&-Q&0\\Q&{\frac {6}{5}}N&{\frac {1}{10}}Nl&-Q&-{\frac {6}{5}}N&{\frac {1}{10}}Nl\\0&{\frac {1}{10}}Nl&{\frac {2}{15}}{Nl}^{2}&0&-{\frac {1}{10}}Nl&-{\frac {1}{30}}{Nl}^{2}\\0&-Q&0&0&Q&0\\-Q&-{\frac {6}{5}}N&-{\frac {1}{10}}Nl&Q&{\frac {6}{5}}N&-{\frac {1}{10}}Nl\\0&{\frac {1}{10}}Nl&-{\frac {1}{30}}{Nl}^{2}&0&-{\frac {1}{10}}Nl&{\frac {2}{15}}{Nl}^{2}\end{array}}\right]$

( 40)

donde K^e es la matriz de rigidez tangente, ${\textstyle {\boldsymbol {K}}_{d}^{e}}$ la matriz de rigidez material y ${\textstyle {\boldsymbol {K}}_{g}^{e}}$ la matriz de rigidez geometrica del elemento de viga 2D. Cabe destacar que Turner et al.[12] fueron los primeros autores a obtener la matriz de rigidez material para un elemento de viga de Timoshenko sin bloqueo por cortante y que en Przemieniecki [13] se describe la matriz de rigidez material en función del coeficiente Φ, definido en (28), a través de la solución de una ecuación diferencial ordinária teniendo en cuenta las hipotesis de Timoshenko. Por otro lado, asumiendo Q  = 0 en la matriz de rigidez geometrica en (40b) se llega a la matriz de rigidez geométrica estandar del elemento de viga de Euler-Bernoulli que se describe como

${\boldsymbol {K}}_{g}^{e}={\frac {N}{30l}}\left[{\begin{array}{cccccc}0&0&0&0&0&0\\0&36&3l&0&-36&3l\\0&3l&4l^{2}&0&-3l&-l^{2}\\0&0&0&0&0&0\\0&-36&-3l&0&36&-3l\\0&3l&-l^{2}&0&-3l&4l^{2}\end{array}}\right]$

( 41)

Esta matriz aparece por primera vez en el artículo publicado por Martin [14] en los años 60.

4. Ejemplos numéricos

Para realizar el análisis no lineal geométrico de los ejemplos presentados en este capítulo, se utilizó un programa escrito en lenguaje Fortran90 denominado co_rotating_2Dbeam.f90 escrito por los autores de este artículo. La tabla 1 resume los principales pasos a ser implementados en un programa computacional para el cálculo tanto del vector de fuerzas internas cuanto de la matriz de rigidez tangente del elemento de viga 2D expresados en coordenadas globales. En todos los ejemplos presentados a continuación fue adoptada una tolerancia para convergencia de 10⁻⁵ . El objetivo en este ítem es mostrar el desempeño del elemento de viga unificado desarrollado en este trabajo para reproducir el comportamiento fuertemente no lineal de algunos sistemas estructurales. Es importante destacar, como será visto en los siguientes ejemplos, que el elemento de viga unificado aquí descrito no presenta bloqueo por cortante, además de la habilidad del mismo en tratar grandes rotaciones de cuerpo rígido.

Tabla 1. Algorítmo de la descripción corrotacional

4.1. Viga en voladizo sometida a flexión pura

Este ejemplo ilustra la habilidad de la formulación corrotacional desarrollada en este trabajo en tratar con grandes rotaciones de sólido rígido debido a la aplicación de la función módulo definida en la ecuación (8) . Se trata de una viga empotrada sujeta a un momento flector aplicado en su extremo libre que se irá enrollando en círculos de radios cada vez más pequeños a medida que se aumenta el valor del momento flector aplicado.

El primer círculo se forma cuando l  = 2πr , o sea, cuando se completa una vuelta. Así, el radio de curvatura será r  = l /2π , la curvatura κ  = 2π /l y el momento flector M  = 2πEI /l . En la figura 7 a se detallan las propiedades mecánicas y geométricas de la viga, la misma fue discretizada en 10, 20 y 40 elementos de viga unificado. Fueron utilizadas mallas más refinadas para representar mejor los círculos de menores radios. Fue impuesto en el extremo libre de la viga una rotación de 2π /10, por lo tanto para completar una vuelta fueron necesarios 10 pasos de carga. En este ejemplo fueron simuladas 8 vueltas con 80 pasos de carga, después de 8 vueltas el radio asume el valor de 1000/16π  = 19.89. Es importante destacar que después de 8 vueltas el extremo libre de la viga sufre una rotación de sólido rígido en el valor de 2880^o . La deformación de la viga mostrada en la figura 7 b fue obtenida utilizando 40 elementos, y se puede observar en esta figura que la viga completa 5 vueltas lo que resulta en la formación de 5 círculos. No fue posible representar con precisión 8 círculos con la malla de 40 elementos porque la longitud de cada elemento es 1000/40 = 25 > 19.89. Por lo tanto, para representar 8 círculos, es necesaria una discretización por encima de 50 elementos. En la figura 7 c se muestran los desplazamientos horizontal y vertical de la extremidad libre de la viga para 8 vueltas con las diferentes mallas utilizadas. Se puede notar que a medida que se refina la malla, mejor es la representación de la magnitud de estos desplazamientos con el aumento de las vueltas. Para la discretización con 10 elementos la trayectoria de equilibrio fue obtenida con un número medio de iteraciones de 5.68, mientras que con 20 elementos este valor fue de 5.04 y con 40 elementos el valor encontrado fue 6.

Figura 7. Viga en voladizo sometida a flexión pura. a) Propiedades geométricas y mecánicas. b) Deformada. c) Trayectorias de equilibrio.

4.2. Viga Williams toggle

Este ejemplo consiste en dos vigas planas inclinadas unidas rígidamente y con los extremos empotrados, sometidas a una carga concentrada P en el vértice. En la figura 8 a se detallan las propiedades geométricas y mecánicas de estas vigas. Cada viga fue discretizada con 10 elementos de viga unificado. En la figura 8 b se muestran dos trayectorias de equilibrio donde fue monitoreado el desplazamiento vertical del vértice para dos condiciones de contorno. La primera condición se refiere a los extremos de ambas vigas empotradas mientras que para la segunda condición se consideró apoyos con dos restricciones. Se compararon estos resultados con los obtenidos por Chan y Chui [15] donde se puede observar una buena concordancia. Para obtener las trayectorias de equilibrio para ambas condiciones de contorno fue controlada la componente vertical del desplazamiento del vértice con el valor de -0.025 para 30 pasos de carga. El numero medio de iteraciones para la condición de contorno empotrada-empotrada fue de 2.93 mientras que el otro caso fue de 3.1. Se puede notar que las trayectorias de equilibrio mostradas en la figura 8 b para ambas condiciones de contorno presentan dos limit points .

Figura 8. Viga Williams toggle. a) Propiedades geométricas y mecánicas. b) Trayectorias de equilibrio.

4.3. Arco circular biarticulado rebajado

En este ejemplo se considera un arco rebajado de directriz circular con ambos extremos apoyados sometido a una carga concentrada P con una excentricidad de 0.2 m en relación a su vértice. De manera que, se consideró la carga P y un momento flector M  = 0.2P aplicados en el vértice del arco. En la figura 9 a se detallan las propiedades geométricas y mecánicas del arco que fue discretizado con 20 elementos de viga unificado. En la figura 9 b se muestra la curva de carga P   versus el desplazamiento vertical ${\textstyle v}$ del vértice. En la figura 9 c se muestra la curva de carga P versus el desplazamiento horizontal u del vértice. Fueron comparados los resultados de ambas trayectorias de equilibrio con los obtenidos por Chan y Chui [15] donde se puede observar una buena concordancia, excepto, pequeñas discrepancias en el tamaño de los lazos en relación al eje de las abscisas presentados en la figura 9 c. Esas discrepancias, probablemente, se deben a que estos puntos fueron obtenidos después de la digitalización de los gráficos y al ajuste de las escalas de los ejes coordenados. Para obtener las trayectorias de equilibrio fue utilizado el método de control variable de los desplazamientos presentados en la referencia [16] . El valor de la longitud de arco fue de 0.0085 para 315 pasos de carga. El número medio de iteraciones fue de 3.6. Se puede notar que las trayectorias de equilibrio del arco mostrado en las figuras 9 b y 9 c son fuertemente no lineales con la presencia de algunos limit points , turning points y loopings .

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

4.4. Pórtico de Lee

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 y . Las condiciones de carga, así como las propiedades geométricas y mecánicas están descritas en la figura 10 a. Este pórtico fue discretizado con 10 elementos de viga unificado. En la figura 10 b se muestran las trayectorias de equilibrio que 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. Los resultados de ambas trayectorias de equilibrio fueron comparados con los obtenidos por Fuji et al.[17] y fue observada una buena concordancia. 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 utilizada fue de 15.1 para 200 pasos de carga. Se presentaron 42 cortes automáticos debido a la divergencia del proceso iterativo. El número medio de iteraciones fue de 4.52. Se puede notar que las trayectorias de equilibrio del Pórtico Lee mostradas en la figura 10 b, también son fuertemente no lineales con la presencia de algunos limit points , turning points y loopings .

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

5. Conclusiones

En este trabajo se mostró, a través de los ejemplos numéricos, la habilidad del elemento finito de viga unificado de tratar con grandes rotaciones de cuerpo rígido. Este elemento integra las teorías de Euler-Bernoulli y de Timoshenko y no presenta bloqueo por deformación de cortante, por lo tanto no produce una sobre rigidez artificial en el comportamiento no lineal de diferentes tipologías estructurales. Fue descrita de manera detallada la relación entre la formulación corrotacional y los modos de deformaciones naturales. La gran ventaja del uso de los modos deformacionales naturales es que los mismos definen a priori los esfuerzos internos auto-equilibrados, lo que lleva a obtener el vector de fuerzas nodales y la matriz de rigidez tangente consistente, descritos en coordenadas globales a través de operaciones algebraicas bastante simples. Por otro lado, en este trabajo se describió de manera explícita el movimiento de sólido rígido del elemento de viga unificado y la obtención de la matriz de rigidez corrotacional. 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.

References

1. [1] B. Skallerud, K. Holthe, B. Haugen; Thin shell and surface crack elements for simulation of combined failure modes; Computer Methods in Applied Mechanics and Engineering, 194/21 (2006), pp. 2619–2640
2. [2] L.L. Yaw, N. Sukumar, S.K. Kunnath; Meshfree co-rotational formulation for two-dimensional continua; International Journal for Numerical Methods in Engineering, 79 (8) (2009), pp. 979–1003
3. [3] C.A. Felippa, B. Haugen; A unified formulation of small-strain corotational finite elements: I. Theory; Computer Methods in Applied Mechanics and Engineering, 81 (2005), pp. 131–150
4. [4] M. Mostafa, M.V. Sivaselvan, C.A. Felippa; A solid-shell corotational element based on ANDES, ANS and EAS for geometrically nonlinear structural analysis; International journal for numerical methods in engineering, 95 (2013), pp. 145–180
5. [5] W.T. Matias e, L.M. Bezerra; Uma abordagem unificada da formulaçao co-rotacional para elementos de treliça 2D, treliça 3D e viga 2D; Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, 25 (2) (2009), pp. 163–190
6. [6] T.N. Le, J.M. Battini, M. Hjiaj; Efficient formulation for dynamics of corotational 2D beams; Computational Mechanics, 48 (2) (2011), pp. 153–161
7. [7] T.N. Le, J.M. Battini, M. Hjiaj; A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures; Computer Methods in Applied Mechanics and Engineering, 269 (2014), pp. 538–565
8. [8] J. Argyris, H. Balmer, J.St. Doltsinis, et al.; Finite element method - the natural approach; Computer Methods in Applied Mechanics and Engineering, 17/18 (1979), pp. 1–106
9. [9] J. Argyris, H.O. Hilpert, G.A. Malejannakis, D.W. Scharpf; On the geometrical stiffness of a beam in space - A consistent V.W. approach; Computer Methods in Applied Mechanics and Engineering, 20 (1979), pp. 105–131
10. [10] S. Krenk; A general format for curved and nonhomogeneous beam elements; Computers & Structures, 50 (1994), pp. 449–454
11. [11] S. Krenk; Non-linear modeling and analysis of solids and structures; Cambridge University Press (2009)
12. [12] M.J. Turner, R.W. Clough, H.C. Martin, L.J. Topp; Stiffness and deflection analysis of complex structures; Journal of the Aeronautical Sciences, 23 (1956), pp. 805–824
13. [13] J.S. Przemieniecki; Theory of matrix structural analysis; McGraw-Hill (1968)
14. [14] H.C. Martin; On the derivation of stiffness matrices for the analysis of large deflection and stability problems, AFFDL-TR-66-80; Air Force Institute of Technology (1966), pp. 697–716
15. [15] S.L. Chan, P.P.T. Chui; Non-linear static and cyclic analysis of steel frames with semi-rigid connections; Elsevier, Amsterdan (2000)
16. [16] 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) (2002), pp. 549–572
17. [17] F. Fujii, K.K. Choong, S.X. Gong; Variable displacement control to overcome turnning points of nonlinear elastic frames; Computers & Structures, 44 (1/2) (1992), pp. 133–136