Pág.1
Elemento Finito Viga
4.
Flexión de Vigas Esbeltas
En este capítulo abordaremos el fenómeno de flexión que experimentan elementos estructurales denominados vigas, los cuales estudiaremos basándonos en las hipótesis de Euler - Bernoulli, el cual desprecia el efecto del esfuerzo cortante generado por la
deformación de corte de la Sección en estudio. Así mismo, veremos algunos aspectos interesantes, cuando estudiemos los elementos finitos tipo viga, como las condiciones de continuidad de las funciones de forma que tienen que satisfacer una continuidad de clase C 1, y los puntos óptimos para el cálculo de los esfuerzos.
4.1. Teoría de Euler -Bernoulli 4.1.1 Teoría básica Consideremos una viga, de longitud L , sección transversal de área A y módulo de inercia I . Sobre dicha viga actúan fuerzas distribuidas, fuerzas puntuales, así como momentos distribuidos y momentos puntuales en un plano xz (Fig. 1.1)
F Mi
q(x)
m(x)
b
y, v
h
x, u L
z , w
z , w
Fig. 4.1.1 Viga de Euler-Bernoulli
La teoría de vigas clásicas de Euler- Bernoulli, se basa en tres hipótesis, además despreciaremos para este estudio la acción de los momento distribuidos m ( x ) = 0 , las
hipótesis son las siguientes: 1) Los desplazamientos verticales de todos los puntos de una sección transversal son pequeños e iguales iguales a los del eje neutro neutro de la viga. 2) El desplazamiento lateral (según eje y de la Fig. 1.1) es nulo. 3) Las secciones transversales normales al eje de la viga permanecen planas antes y después de la deformación.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
1
Pág.2
Elemento Finito Viga
F Mi
q(x) A
,
Z
B ,
x, u
w( x) A B
,
,
θ
=
,
A
dw dx
θ B
,
B
z, w
u ( x)
Fig. 4.1.2 Desplazamientos de la sección de una viga
Según las hipótesis establecidas el campo de desplazamiento de un punto cualquiera se puede escribir escribir como:
u( x, y , z ) = − z
dw dx
v ( x, y , z ) = 0
w( x, y , z ) = w( x )
(1.1) (1. 1)
La deformación en un punto material se obtiene por las respectivas definiciones:
du
ε x
=
ε y
= ε z = γ xy = γ xz = γ yz = 0
dx
El único esfuerzo no nulo
= − z
d 2 w
σ x se
dx 2
relaciona con su correspondiente deformación
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.2)
ε x por:
2
Pág.3
Elemento Finito Viga
σ x
= E ε x = − zE
2
d w dx 2
(1.3)
Se define el momento flector positivo M (según Fig.1.3) de una sección como:
M =
− ∫∫ zσ x dA = ∫∫ z E 2
A
A
d 2 w dx
2
dA = EI
d 2 w dx
2
= EI χ
(1.4)
Donde I es el momento de inercia de la sección transversal con respecto al eje y siendo χ la
M
curvatura del eje de la viga.
σ x
σ x M x
− σ x
− σ x
z
Fig. 4.1.3 Convenio de signos para los esfuerzos y el momento flector M .
El equilibrio de la viga está dado por el principio de los trabajos virtuales ( PTV ), que se escribe de la siguiente forma:
∫
L
∫
δ ε xσ x dV = δ wq( x ) dx + 0
V
⎞ ⎟ ∑ δ w F + ∑δ ⎛ ⎜⎝ dw dx ⎠ q
P
M j
1 i
i =1
j =1
(1.5)
j
La integral sobre el volumen de la viga del primer miembro representa el trabajo generado por la deformación virtual y puede simplificarse de la siguiente forma:
⎡ ⎤ ⎛ d 2 w ⎞ 2 ∫V δ ε xσ x dV = ∫0 ⎢⎣ A∫ − z E ε x dA⎥⎦δ ⎜⎜⎝ dx 2 ⎠⎟⎟dx = L
L ⎛ d 2 w ⎞ ⎛ d 2 w ⎞ = ∫ δ ⎜⎜ 2 ⎟⎟ EI ⎜⎜ 2 ⎟⎟ dx = ∫ δ ( χ ) Mdx dx ⎠ ⎝ dx ⎠ 0 ⎝ 0 L
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.6)
3
Pág.4
Elemento Finito Viga
El trabajo de deformación virtual, es la integral sobre la longitud de la viga del producto del momento flector por la correspondiente curvatura virtual.
4.1.2 Discretización en elementos finitos de dos nodos La ecuación de equilibrio dada por el PTV , indica que solamente es necesario determinar el campo de desplazamientos verticales w( x ) , para satisfacer el equilibrio. Sin embargo se observa que el trabajo virtual interno depende de la segunda derivada de la flecha w (curvatura
χ ),
lo cual implica que se tiene que utilizar elementos con
continuidad de clase C 1 (la variable y su primera derivada han de ser continuas) para evitar singularidades en el cálculo de las integrales. Esta condición puede interpretarse físicamente, pues la primera derivada
dw dx
coincide con la pendiente de la deformada de
la viga. Por tanto, dicha derivada debe ser continua para garantizar que la deformada sea una curva suave. El elemento más sencillo de clase C 1 es el unidimensional de dos nodos. La condición de continuidad de la primera derivada de la flecha, obliga a tomar dw
el giro
dx
como variable. Por consiguiente el número total de variables nodales del
dw ⎞ elemento es cuatro ( wi y ⎛ ⎜ ⎟ por nodo)
⎝ dx ⎠ i
(
w1
w2
1
2
dw dx
)1
z
L
(
x dw dx
(e)
)2
(e)
Fig. 4.1.4 Variables en el elemento de dos nodos
Las variables anteriores definen perfectamente una variación cúbica de la flecha w( x ) = a0
+ a1 x + a2 x 2 + a3 x 3
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.7) 4
Pág.5
Elemento Finito Viga
Las constantes
ai
se calculan sustituyendo adecuadamente los valores de la flecha y sus
derivados en los nodos, lo que genera un sistema de ecuaciones algebraicas con cuatro incógnitas.
w1
= a0 + a1 x1 + a2 x12 + a3 x13
⎛ dw ⎞ = a + 2a x + 3a x 2 ⎜ ⎟ 1 2 1 3 1 ⎝ dx ⎠1 w2
= a0 + a1 x2 + a2 x22 + a3 x23
⎛ dw ⎞ = a + 2a x + 3a x 2 ⎜ ⎟ 1 2 1 3 2 ⎝ dx ⎠ 2
(1.8)
Una vez resuelto el sistema anterior se puede escribir como: L(
e)
e dw ⎞ L( ) ⎛ dw ⎞ ⎛ w(ξ ) = N 1 (ξ ) w1 + N 1 (ξ ) ⎜ ⎟ + N 2 (ξ ) w2 + N 2 ⎜ ⎟ (1.9) 2 ⎝ dx ⎠1 2 ⎝ dx ⎠ 2
Para lo cual se utilizan los funciones de forma hermíticas en un sistema de coordenadas parametrizado ξ ∈ [− 1, 1] .
ξ
ξ = −1
ξ = +1
L
Fig. 4.1.5 Coordenadas parametrizadas
En este sistema parametrizado la flecha se aproxima como:
⎛ dw ⎞ ⎛ dw ⎞ ⎟⎟ + N 2 (ξ )w2 + N 2 (ξ ) ⎜⎜ ⎟⎟ ⎝ d ξ ⎠1 ⎝ d ξ ⎠ 2
w (ξ ) = N 1 (ξ )w1 + N 1 (ξ ) ⎜⎜
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.10)
5
Pág.6
Elemento Finito Viga
N 1
N 1
1
1
N 2
1
1 N 2
Fig. 4.1.6 Funciones Hermíticas
Las funciones de forma hermíticas están dadas por:
Con
N 1
=
1 (2 − 3ξ + ξ 3 ) ; N 2 = 1 (2 + 3ξ − ξ 2 ) 4 4
N 1
=
1 (2 − 3ξ + ξ 3 ) ; N 2 = 1 (− 1 − ξ + ξ 2 + ξ 3 ) 4 4
ξ =
2
( x − x m ) y ( e)
L
xm
=
x1 + x2
2
(1.11)
(1.12)
La ecuación (1.9) puede escribirse haciendo uso del cálculo matricial de la siguiente manera: w=
Nu ( e )
(1.13) T
donde:
⎡ ⎡ ⎛ dw ⎞ L( e ) L( e ) ⎤ ⎛ dw ⎞ ⎤ (1.14) (e) N = ⎢ N 1 , N 1 ; N 2 ; N 2 u y = ; ; ; w w ⎜ ⎟ ⎟ ⎥ 2 ⎜ ⎥ ⎢ 1 dx 2 2 dx ⎝ ⎠ ⎝ ⎠ 1 2⎦ ⎣ ⎦ ⎣
Son las matrices de las funciones de forma y el vector de desplazamientos nodales del elemento. Las funciones de forma están aproximados por polinomios hermíticos la figura 1.6 muestra las funciones de forma hermíticas , se puede observar que las funciones N 1 y
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
6
Pág.7
Elemento Finito Viga
N 2 valen la unidad en un nodo y cero en el otro, mientras que los giros son cero para
ambos sucediendo lo contrario con las funciones N 1 y N 2 De la ecuación (1.12) se deduce que
dx =
( )
L e
2
d ξ ;
dx d ξ
dw dx
=
=
( e)
L
, con lo que:
2
2 dw (e)
L
d ξ
2
d w
y
dx
2
=
4 d 2 w ( e )2
L
dq
2
(1.15)
Con la ecuación (1.9) y (1.15) se pueden determinar la pendiente y curvatura en un punto del elemento de coordenada ξ .
χ =
2
d w dx
2
(e) ⎛ d 2 N 1 L = ( e ) 2 ⎜⎜ 2 w1 + 2 ( L ) ⎝ d ξ
4
2
d N 1 ⎛ dw ⎞ d ξ
⎜ ⎟ + ⎝ dx ⎠1
2
d w 2
d ξ
w2
+
( e)
L
2
2
d N 2 d ξ
⎛ dw ⎞ ⎞⎟ = ⎜ ⎟ ⎟ ⎝ dx ⎠ 2 ⎠
⎧ w1 ⎫ ⎪ ⎛ dw ⎞ ⎪ ⎡ − 6ξ (3ξ − 1) − 6ξ (3ξ + 1) ⎤ ⎪⎪ ⎜⎝ dx ⎠⎟1 ⎪⎪ = ⎢ ( e ) 2 , ( e ) , ( e) 2 , ( e ) ⎥ ⎨ = B f u ( e) ⎬ L ( L ) L ⎦⎥ ⎪ w2 ⎪ ⎣⎢ ( L ) dw ⎞ ⎪ ⎪⎛ ⎜ ⎪⎩⎝ dx ⎠⎟ 2 ⎪⎭
(1.16)
Siendo B f la matriz de deformación de flexión del elemento. Finalmente la expresión de los trabajos virtuales de un elemento toma la siguiente forma:
∫
L( e )
⎛ +1 ( e ) T T L ⎞ δχ EI χ dx = ⎜ ∫ [δ u ] B f ( EI )B f d ξ ⎟ u ( e ) = ⎜ 2 ⎠⎟ ⎝ −1 +1
=
∫ [δu ] N
−1
( e ) T
( e)
T
q( x )
L
2
d ξ +
2 sw ⎞ δ + δ⎛ w F ⎜ ⎟ M j ∑ ∑ i i i =1 j =1 ⎝ dx ⎠ j 2
(1.17)
operando se obtiene en forma matricial:
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
7
Pág.8
Elemento Finito Viga
K ( e ) u ( e ) − q ( e ) = f ( e)
(1.18)
donde la matriz de rigidez se calcula en forma explícita por:
K (e)
⎡ 12 6 L( e ) 6 L( e ) ⎤ − 12 (e) ⎢ +1 ( e) 2 ( e) ( e) 2 ⎥ (e) 4 6 2 ( ) − ( ) ⎥ L L L L EI ⎛ ⎞ = ∫ B T f B f EI d ξ = ⎜ 3 ⎟ ⎢ 2 12 − 6 L( e ) ⎥ ⎝ L ⎠ ⎢ −1 ⎢ 2⎥ 4( L( e ) ) ⎦ ⎣ sim
(1.19)
La matriz de rigidez coincide con las obtenidas haciendo uso de las clásicas ecuaciones de resistencia en el cálculo matricial de estructuras. Esto se debe a que las funciones hermíticas que aproximan al desplazamiento en el elemento de dos nodos, coincide
exactamente con la que se obtiene integrando la ecuación diferencial de equilibrio de la viga. Por otra parte, el vector de fuerzas nodales equivalentes debido a una fuerza uniformemente distribuida de intensidad q sobre el elemento es: +1
q
(e)
= ∫N −1
(e)
T
q
L
2
d ξ = qL
(e)
⎡ 1 L( e ) 1 L( e ) ⎤ ⎢ 2 ; 12 , 2 − 12 ⎥ ⎣ ⎦
(1.20)
Y el vector de fuerzas nodales de equilibrio f (e) necesario para el ensamblaje
f ( e ) = [V 1 , M 1 , V 2 , M 2 ]T
(1.21)
una vez determinados los desplazamientos y giros nodales, se puede calcular el momento flector en cualquier punto del elemento a través de la siguiente expresión: M = EI χ = EI B f u ( e )
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.22)
8
Pág.9
Elemento Finito Viga
Ejemplo: Calcular los desplazamientos en los extremos y las reacciones en el extremo empotrado de una viga en voladizo bajo la acción de una carga uniformemente distribuida, utilizando un solo elemento.
q EI L
olución por el Método de los Elementos Finitos ( Diagrama de Cuerpo Libre)
q x , u
M1
V1
z , w
Fig. 4.1.7 Viga en voladizo bajo carga uniforme
⎧ qL − V ⎫ 1 ⎪ ⎪ 2 − 12 6 L ⎤ ⎧W 1 ⎫ ⎪ qL2 ⎪ ⎪ ⎪ − 6 L 2 L2 ⎥⎥ ⎪ θ1 ⎪ ⎪⎪ 2 − M 1 ⎪⎪ ⎨ ⎬=⎨ ⎬ 12 − 6 L⎥ ⎪W 2 ⎪ ⎪ qL ⎪ ⎥ 2 ⎪ 4 L2 ⎦ ⎪⎩ θ2 ⎪⎭ ⎪ 2 qL ⎪ − ⎪ ⎪⎩ 12 ⎪⎭
⎡12 6 L ⎢ 4 L2 ⎛ EI ⎞ ⎢ ⎜ 3⎟ ⎝ L ⎠ ⎢ ⎢ ⎣
Sabiendo que w1 y
θ 1 =
(1.23)
0 , por ser un empotramiento, se pueden prescindir de las dos
primeras filas y columnas correspondientes a dichos grados de libertad, obteniendo luego los desplazamientos ( w2 y w2
=
ql
4
8 EI
; θ 2
θ 2 )
del nodo dos.
3 dw ⎞ ql ⎛ =⎜ ⎟ = ⎝ dx ⎠ 2 6 EI
(1.24)
Con estos valores se determinan las reacciones en el empotramiento. Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
9
Pág.10
Elemento Finito Viga
V 1
= ql y
M 1
=
ql
2
2
(1.25)
Dichos valores coinciden con los valores exactos, sin embargo dicha coincidencia ocurre sólo en los nodos ya que la variación exacta de la flecha es un polinomio de cuarto grado, mientras que la utilizada en la aproximación es una de tercer grado. En la figura a continuación se muestra la ley de momentos flectores obtenida por el método de los elementos finitos y la exacta.
χ = EI B f u M MEF = EI M Ex
=
qL2
8
=
ql
6
2
−
ql
2
4
ξ
(1 − ξ )2
(1.26)
ξ
ξ = -1
ξ = +1 ξ = ξ = -
+
1 3
1 3
2
qL
2
Solución MEF (1 Elemento)
Solución Exacta
Fig. 4.1.8 Ley de momentos
Es interesante advertir que amas soluciones coinciden con los dos puntos de Gauss para ξ = ±
1 3
. Esta coincidencia no es fortuita y dada su importancia en la integración
numérica estudiaremos con más detalle en el apartado siguiente.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
10
Pág.11
Elemento Finito Viga
4.1.3 Puntos óptimos para el cálculo de esfuerzos y deformaciones Al obtener los esfuerzos a partir de las derivadas de los desplazamientos su aproximación es siempre de menor orden que la de los desplazamientos. Por lo general, si las funciones de forma son polinomios de grado p, la aproximación de las tensiones será de polinomios de grado p-1 ó p-2, según se obtenga en función de la primera o segunda derivada del campo de desplazamientos. Está demostrado que las tensiones obtenidas en el método de los elementos finitos pueden considerarse como un ajuste por mínimos cuadrados ponderados de la solución exacta de los esfuerzos ( Zienkiewcz). Esta propiedad cobra una vital importancia en la evalución de los esfuerzos dentro de los elementos finitos, que justamente coinciden en los puntos de intersección con la solución exacta . La dificultad surge
que en la mayoría de los casos no se conoce una distribución real del campo de esfuerzos. Esta dificultad puede sortearse haciendo uso de la siguiente propiedad de la integración numérica de Gauss – Legendre: en los puntos de una cuadratura de Gauss- Legendre de orden n, un polinomio de grado n y otro de grado n-1, obtenido como ajuste por mínimo cuadrado del anterior toman el mismo valor.
Aclaremos esta propiedad en dos
ejemplos.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
11
Pág.12
Elemento Finito Viga
Ejemplo 1: comprobar que un polinomio de segundo orden y otro de primer orden obtenido por ajuste por mínimo cuadrado del anterior se cortan en los puntos de la cuadratura de Gauss-Legendre de orden 2. Tomemos un polinomio de segundo grado ( n = 2) f ( x ) = 1 + x + x 2
y obtengamos el de primer grado ( n = 1) que ajusta a
f ( x )
por mínimos cuadrados. Tal
que: g ( x ) = a + bx
y de forma que: +1
M =
∫ [ f ( x ) − g ( x)]
2
−1
+1
dx
= ∫ [(1 − a ) + (1 − b) x + x 2 ] dx 2
−1
Sea un mínimo Las constantes a y b se obtienen de las siguientes condiciones:
∂ M =0 ⇒ ∂a
+1
∂ M =0 ⇒ ∂b
+1
y resolviendo se obtiene
∫ − 2[(1 − a ) + (1 − b) x + x ] dx = 0 2
−1
∫ − 2b[(1 − a ) + (1 − b) x + x ] dx = 0 2
−1
a
=
4 3
y b =1.
En la figura 1.9, se muestran las funciones f ( x) y g( x). Puede observarse que ambos polinomios se cortan en los puntos de la cuadratura de Gauss-Legendre de orden 2.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
12
Pág.13
Elemento Finito Viga
y
f(x) = 1 + x + x 2
g(x) =
-1
1 x = 3
1 x = + 3
4 + x 3
x
+1
Fig. 4.1.9 Puntos de intersección de los polinomios f ( x ) y g ( x )
Ejemplo 2: comprobar que un polinomio cúbico y otro cuadrático obtenido por ajuste por mínimos cuadrados del anterior se cortan en los puntos de la cuadratura de Gauss Legendre de orden 3. (se deja al lector la solución de este ejemplo, los puntos de la
cuadraturas de Gauss-Legendre son: x1 = -0.7746; x2 = 0; x3 = 0.7746) De lo anterior pueden realizarse las dos siguientes conclusiones: 1) Si la distribución exacta de los esfuerzos
σ
(o las deformaciones
ε ),
es un
polinomio de grado n, y la aproximada obtenida por elementos finitos es de grado n-1, la evaluación de los esfuerzos σ o de las deformaciones ε en los puntos de la
cuadratura de Gauss-Legendre de orden n es exacta. 2) Si los polinomios que representan las soluciones exactas y de elementos finitos para σ
y
ε difieren
en más de un grado, la evaluación de
σ
y
ε en
los puntos de la
cuadratura de Gauss-Legendre aproxima un término más del desarrollo de la serie de Taylor de la solución exacta que en cualquier otro punto del elemento. Resumiendo se puede decir que los puntos de la cuadratura de Gauss-Legendre tienen la interesante propiedad de aproximar con un orden mayor los esfuerzos (y las deformaciones) y por consiguiente deben evaluarse en dichos puntos y a partir de los valores obtenidos proceder si se desea, a la extrapolación a los nodos. Por ello dichos puntos de consideran óptimos para el cálculo de los esfuerzos y deformaciones. Los
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
13
Pág.14
Elemento Finito Viga
conceptos anteriormente mencionados son rigurosamente válidos para elementos unidimensionales. Asimismo, se ha comprobado que la utilización de las mismas ideas para la integración numérica de Gauss-Legendre en dos y tres dimensiones lleva a una mejora sustancial en la evaluación de los esfuerzos y deformaciones. En la Fig. 1.10 se muestran los puntos óptimos para el cálculo de esfuerzos y deformaciones en algunos de los elementos uni y bidimensionales.
Barra a tracción
Flexión de vigas
C 1
x
x
x
C 0 {
x
x
x x
x
Elasticidad 2D
X
x
x
X
x
X
x
x
x
x O
x
x
x
x
Placas y Láminas
C 0
{
x x
x x
X
x
x
C 1
x
x
x
x
x
x x
x
x O
x
x
x
Fig. 4.1.10 Puntos óptimos para la evaluación de los esfuerzos
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
14
Pág.15
Elemento Finito Viga
4.2 FLEXIÓN DE VIGAS DE TIMOSHENKO 4.2.1 Teoría básica La nueva hipótesis introducida por Timoshenko es la de considerar que “las secciones planas normales al eje de la viga antes de la deformación, permanecen planos pero no necesariamente normales al eje después de la deformación”.
Fi q(x) x, u
β
θ φ
w( x) dw dx
= β
z , w
Fig. 4.2.1 Teoría de flexión de vigas de Timoshenko
Esta hipótesis represente una mejor aproximación a la deformación real de la sección transversal. Esta hipótesis tiene su sustento para las vigas de gran canto donde las secciones transversales dejan de conservarse después de la deformación. Esta hipótesis asume un giro medio para la sección, de manera que para efectos prácticos se considera plana. De la figura 1.11 se deduce que el giro de la sección se expresa como θ =
donde
dw dx
d w d x
+ φ
(1.27)
es la pendiente de la deformada del eje de la viga y φ un giro adicional
debido a la deformación por cortante. El campo de desplazamiento se expresa como u ( x, y, z ) = − z θ ( x ) v ( x, y, z ) = 0 w( x, y, z ) = w ( x )
(1.28)
De dicho campo de desplazamientos podemos ahora determinar las deformaciones no nulas
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
15
Pág.16
Elemento Finito Viga
ε x
=
d u
γ xz
=
d w
d x
= − z
d x
+
d θ d x
d u d z
=
d w d x
− θ = −φ
(1.29)
Por lo tanto la teoría de Timoshenko equivale a considerar el efecto de la deformación por cortante transversal coincidiendo la magnitud de dicha deformación con el giro adicional de la normal φ . Las deformaciones no nulas generan tensiones que están relacionadas por las respectivas ecuaciones constitutivas: σ x = E ε x
τ xz
d θ
= − zE
d x
= − zE χ
⎛ d w ⎞ = Gγ xz = G⎜⎜ − θ ⎟⎟ d x ⎝ ⎠ d θ
Aquí G es el módulo de corte y χ =
d x
(1.30)
la curvatura del eje de la viga.
El momento flector y el esfuerzo cortante que generan dichas tensiones se defines como
∫
∫
M = − z σ x d A = z E A
A
2
d θ d x
χ dA = EI
con dA = d y d z
∫
⎛ d w ⎞ − θ ⎟⎟ = GAγ xz d x ⎝ ⎠
Q = τ xz d A = GA ⎜⎜ A
(1.31)
σ x Q
M x
z Tensión Normal σ x Distribución supuesto=Distri. exacta
Q x
τ xz
z Tensión tangencial τ xz Distri. supuesto
x
τ xz
z Tensión Tangencial τ xz Distri. exacta
Fig. 4.2.2 Teoría de vigas de Timoshenko. Distribución de tensiones normales y tangenciales
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
16
Pág.17
Elemento Finito Viga
Obsérvese que la variación de σ x en el canto es lineal, lo cual puede considerarse como “exacto ” dentro de la hipótesis de la teoría de vigas. Por el contrario, la variación de la tensión tangencial τ xz en el canto se supone constante, lo cual está en clara contradicción con la distribución polinómica de la teoría de vigas. Para sortear este problema se acepta la hipótesis de tensión tangencial constante pero modificada por un coeficiente de manera que el trabajo de deformación de tensión tangencial constante coincida con el “exacto” de la teoría de vigas. Así se toma
= α Gγ xz
τ xy
y Q = α AGγ xz
= A ∗ Gγ xz
(1.32)
donde α es el coeficiente de forma o distorsión de la sección y A ∗ = α A se denomina área reducida. Con las tensiones, deformaciones y desplazamientos antes planteados podemos ahora expresar el principio de los trabajos virtuales ( PTV ) de la siguiente forma L
∫ (δ ε σ x
x
p
+ δγ xzτ xz ) dV = ∫ δ wqd x + ∑δ wi F i + ∑δ θ j M j i =1
0
v
q
(1.33)
j =1
El primer miembro de la ecuación anterior puede modificarse como:
⎡ ⎛ d θ ⎞ ⎛ d w ⎞⎤ ⎜ ⎟ − + σ δ τ δ z ∫V ⎢⎣ x ⎜⎝ d x ⎠⎟ xz ⎜⎜⎝ d x − θ ⎠⎟⎟⎥⎦ dV = L
⎡
∫ ⎣⎢
⎛
⎞
⎛
⎞⎤
⎝ A
⎠
⎝ A
⎠⎦⎥
∫
∫
= ⎢δχ ⎜⎜ − z σ x dA ⎟⎟ + δ γ xz ⎜⎜ τ xz dA ⎟⎟⎥ dx = 0
L
∫
= [δχ M + δ γ xz Q ] dx = 0
⎡ ⎛ d θ ⎞ d θ ⎛ d w ⎞ ∗ ⎛ d w ⎞⎤ ⎟ δ + EI ∫0 ⎣ ⎝ d x ⎠⎟ d x ⎜⎜⎝ d x − θ ⎠⎟⎟ GA ⎜⎜⎝ d x − θ ⎠⎟⎟⎥⎦ dx
L
= ⎢δ ⎜⎜
(1.34)
En la ecuación (1.34) se aprecia que en el integrando aparecen únicamente derivadas primeras de las flechas y el giro. Esto exige únicamente su continuidad para garantizar la integrabilidad, lo que permite la utilización de elementos finitos de clase C 0.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
17
Pág.18
Elemento Finito Viga
4.2.2 Elementos Finitos para la flexión de vigas de Timoshenko Primeramente consideraremos el elemento viga de Timoshenko más sencillo de dos nodos. A diferencia de la teoría de vigas de Euler – Bernoulli, la flecha y el giro son variables independientes y de continuidad C 0. Así podemos interpolar por separado cada una de ellas por w (ξ ) = N 1 (ξ ) w1
+ N 2 (ξ )w2
θ (ξ ) = N 1 (ξ )θ 1
+ N 2 (ξ ) θ 2
(1.35)
1 1 ;y N 2 = (1 + ξ ) y además w1 , θ 1 y w2 , θ 2 son las 2 2 flechas y giros de los nodos 1 y 2 del elemento. con
N 1
= (1 − ξ )
Haciendo uso de la ecuación (1.35) la curvatura se obtiene χ =
d θ d x
=
d ξ d θ
⋅
d x d ξ
d ξ ⎡ d N 1
=
⎢ d x ⎣
d ξ
θ 1
dN 2
+
d ξ
⎤ ⎦
θ 2 ⎥
(1.36)
y la deformación de cortante (o cizalladura) γ xz
=
d w d x
− θ =
d ξ ⎡ dN 1
⎢ d x ⎣
d ξ
w1
+
⎤ ⎦
dN 2
w2 ⎥ − ( N 1θ 1
d ξ
+ N 2θ 2 ) d ξ
Utilizando una formulación isoparamétrica se tiene que
d x
=
(1.37)
2 L(
e)
y con ello las
ecuaciones (1.36) y (1.37) pueden ser escritas en forma matricial como χ = B f u (e ) γ xz
= B c u (e )
(1.38)
donde
⎡ 2 B f = ⎢0, (e ) ⎣ L ⎡ 2 B c = ⎢ (e ) ⎣ L
dN 1 d ξ
dN 1 d ξ
; 0;
2 dN 2 ⎤ (e )
L
; − N 1 ;
1 1 = ⎡⎢0; − ( e ) ; 0 ; ( e) ⎤⎥ ⎥ d ξ ⎦ ⎣ L L ⎦
2 dN 2 (e)
L
d ξ
⎤ ⎡ ⎦ ⎣
; − N 2 ⎥ = ⎢ −
1 (e)
L
;−
(1 − ξ ) 1 2
;
;− (e)
L
(1 + ξ ) ⎤ 2 ⎥⎦
(1.39)
Son las matrices de deformación de flexión y cortante del elemento, y
u (e ) = [w1 ; θ 1 ; w2 ; θ 2 ]T
(1.40)
es el vector de desplazamientos nodales.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
18
Pág.19
Elemento Finito Viga
La expresión de los trabajos virtuales (1.33) puede escribirse como
⎛
[δ u ( ) ] ⎜⎜ ∫ [B ( EI ) B e
T
⎝ L( )
T f
⎞ + B T c (GA∗ )B c ]dx ⎟⎟ u (e ) = ⎠
f
e
⎛
⎞
= [δ u (e ) ] ⎜⎜ N T (q ) d x + f (e ) ⎟⎟ T
∫ ⎝
(1.41)
⎠
L( e )
y tras simplificar los desplazamientos virtuales queda
K f (e ) + K (ce ) u (e ) − q ( e ) = f ( e )
(1.42)
donde
K (e ) = K f (e ) + K (ce ) y
K f (e ) = ∫ B T f ( EI ) B f d x ; K (ce ) = ∫ B T c (GA∗ ) B c d x (e )
(1.43)
(e )
L
L
son las matrices de rigidez correspondientes a los efectos de flexión y cortante cuya suma es la matriz de rigidez total del elemento;
q (e ) = ∫ N T qd x ;
con
N = [ N 1 ; 0 ; N 2 0]
(1.44)
L( e )
es el vector de fuerzas nodales equivalentes debidas a la carga distribuida q; y
f (e ) = [V 1 ; M 1 ; V 2 ; M 2 ]T
(1.45)
el vector de fuerzas nodales de equilibrio que permite ensamblar las contribuciones de los distintos elementos en la matriz de rigidez y el vector de fuerzas globales. Todas las integrales anteriores pueden transformarse sobre el dominio normalizado del elemento. Así teniendo en cuenta que d x = 1
K f = ∫ B ( EI ) B f (e )
T f
−1
L(
e)
2
d ξ ;
L(e )
2
d ξ se tiene que: 1
Kc = ∫B (e )
−1
T c
(GA ) B
L(
∗
c
e)
2
d ξ
(1.46)
y 1
q = ∫N (e )
−1
(e )
L q
T
2
d ξ
(1.47)
Expresiones que pueden evaluarse numéricamente por una cuadratura unidimensional de Gauss-Legendre.
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
19
Pág.20
Elemento Finito Viga
4.2.3 Efecto del bloqueo de la solución
Para determinar los valores exactos de la matriz de rigidez a la flexión K f es necesario un solo punto de integración ya que todos los términos del integrando son constantes, obteniendo finalmente
⎡0 0 (e ) ⎢ EI ⎞ ⎢0 1 ⎛ (e ) K f = ⎜ ⎟ ⎝ L ⎠ ⎢0 0 ⎢ ⎣0 − 1
0 0⎤ 0 − 1⎥⎥ 0 0⎥ ⎥ 0 1⎦
(1.48)
Para determinar la matriz de rigidez de cortante es necesario la utilización de dos puntos de integración para hallar los valores exactos por aparecer en el integrando de K c términos de segundo orden en ξ (debido a los productos N i N j ), obteniéndose
(e )
Kc
⎡ ⎢1 ⎢ (e ) ⎢ ⎛ GA∗ ⎞ ⎢ ⎟⎟ ⎢ = ⎜⎜ L ⎝ ⎠ ⎢ ⎢ ⎢ ⎢ ⎣
(e )
L
−1
2
( L( ) )2 e
3
−
(e )
L
2 1
⎤ 2 ⎥⎥ ( L( e) ) 2 ⎥ 6 ⎥ (e ) L ⎥ ⎥ − 2 ⎥ ( L(e) )2 ⎥ ⎥ 3 ⎦ (e )
L
(1.49)
Para apreciar el efecto de la integración numérica estudiaremos la flexión de la viga en voladizo bajo la acción de una carga puntual. La ecuación matricial del equilibrio global es la siguiente: P A -A
1
A
2
A L
h b
Fig. 4.2.3 Viga en voladizo con un elemento de viga de Timoshenko de 2 nodos
K f (1) + K (c1) u (1) = f
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.50)
20
Pág.21
Elemento Finito Viga
⎡ GA∗ ⎢ L ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
∗
GA ∗
−
2
∗
∗
⎤ ⎥ 2 ⎥ ⎛ GA∗ EI ⎞⎥ ⎧ w1 ⎫ ⎧ V 1 ⎫ ⎜⎜ ⎟⎟ ⎪ ⎪ ⎪ ⎪ L − 6 L ⎝ ⎠⎥ ⎪ θ 1 ⎪ = ⎪ M 1 ⎪ ⎥ ⎨w ⎬ ⎨ P ⎬ ∗ GA ⎥ ⎪ 2⎪ ⎪ ⎪ − 2 ⎥ ⎪⎩θ 2 ⎪⎭ ⎪⎩ 0 ⎪⎭ ∗ ⎛ GA EI ⎞⎥ ⎜⎜ ⎟⎟⎥ L + 3 L ⎝ ⎠⎦
GA
GA
L GA∗
⎛ GA EI ⎞ ⎜⎜ ⎟⎟ − L + 3 2 L ⎝ ⎠
∗
GA L Simétrico
(1.51)
sabiendo que por el empotramiento se tiene la siguientes condiciones de borde w1 = 0 y θ 1 = 0 Eliminando fila y columna correspondiente a los desplazamientos anteriores se obtiene el sistema de ecuaciones simplificado:
⎡ GA∗ ⎢ L ⎢ ∗ ⎢− GA ⎢⎣ 2
∗
⎤ ⎥ ⎧w 2 ⎫ ⎧ P ⎫ 2 ⎥ ⎬=⎨ ⎬ ⎛ GA∗ EI ⎞⎥ ⎨ θ ⎜⎜ ⎟⎟ ⎩ 2 ⎭ ⎩ 0 ⎭ L + L ⎠⎥⎦ ⎝ 3 −
GA
(1.52)
La solución se encuentra por
⎡⎛ L L3 ⎞ ⎜ ⎟ + ⎧w2 ⎫ β ⎢⎜⎝ GA∗ 3 EI ⎠⎟ ⎢ ⎨θ ⎬ = F f = 1 β + L2 ⎢ ⎩ 2⎭ ⎢⎣ EI
⎤ ⎥ EI ⎥ ⎧ P ⎫ ⎨ ⎬ L ⎥ ⎩ 0 ⎭ EI ⎥⎦ L2
donde F = K −1 es la matriz de flexibilidad y β =
12 EI ∗ 2
GA L
(1.53)
. De (1.53) se deduce la flecha
en el extremo libre y vale
⎛ L L3 ⎞ ⎜ ⎟P + w2 = β + 1 ⎜⎝ GA∗ 3 EI ⎠⎟ β
En el caso de una sección rectangular I =
(1.54)
bh 3
12
; A ∗ =
5 bh y ν = 0,25 6
2
⎛ h ⎞ 3 β = 3 ⎜ ⎟ = 2 ⎝ L ⎠ λ
(1.55)
⎛ L ⎞ ⎝ h ⎠
donde λ = ⎜ ⎟ se denomina coeficiente de esbeltez de la viga. La expresión “exacta” de la matriz de flexibilidad de una viga con y sin la inclusión del efecto del esfuerzo cortante de acuerdo con la teoría de vigas clásica es:
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
21
Pág.22
Elemento Finito Viga
a) Sin esfuerzo cortante (Euler-Bernoulli)
b) Con esfuerzo cortante (Timoshenko) 3 2 ⎡⎛ L L ⎞ L ⎤ ⎟⎟ ⎢⎜⎜ ∗ + ⎥ 3 2 EI EI GA ⎠ ⎥ F = ⎢⎝ 2 ⎢ L L ⎥ ⎢⎣ ⎥ 2 EI EI ⎦
⎡ L3 L2 ⎤ ⎢ ⎥ 2 EI ⎥ F = ⎢ 3 EI 2 L ⎥ ⎢ L ⎢⎣ 2 EI EI ⎥⎦
(1.56)
Por lo tanto, la flecha “exacta ” en el extremo de la viga es: a)
Sin esfuerzo cortante
(w2 )exacta = f
L3
3 EI
b) Con esfuerzo cortante
(w2 )exacta C
P
⎛ L L3 ⎞ = ⎜⎜ ∗ + ⎟⎟ P 3 EI GA ⎝ ⎠
(1.57)
Para una viga esbelta (valores de λ elevados) el efecto del esfuerzo cortante es despreciable y la solución numérica obtenida debe coincidir con la expresión (a) de la ecuación (1.57). De la ecuación (1.54) y (1.57) se deduce el cociente entre la solución de elementos finitos y la teoría para vigas esbeltas es
ϕ =
w2
(w2 ) f exacta
⎛ L L3 ⎞ ⎜⎜ ∗ + ⎟P 3 EI ⎠⎟ 3 (4λ 2 + 3) β ⎝ GA = = 2 2 4λ (λ + 3) β + 1 ⎛ L3 ⎞ ⎜⎜ ⎟⎟ P ⎝ 3 EI ⎠
(1.58)
Lógicamente, el valor de ϕ debería tender a la unidad a medida que la esbeltez de la viga aumenta (mayor λ ). ϕ
4 Solución exacta 1 elemento
3 4λ 2
+3
2 elementos
4λ 2
1 elemento
}
Integración reducida de k c (e)
(1 punto)
de k } Integración reducida 2 con
2
puntos
ϕ = 0,938; λ → ∞
1 3(4λ 2 + 3) 4λ 2 (λ 2 + 3)
ϕ = 0,75; λ → ∞ λ
1
2
3
4
5
Fig. 4.2.4 Viga en voladizo analizado con un elemento de Timoshenko de dos nodos. Variación del cociente entre la solución obtenida para la flecha en el extremo y la exacta de la teoría de Euler-Bernoulli con el coeficiente de esbeltez Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
22
c
(e)
Pág.23
Elemento Finito Viga
La figura 4.2.4 nos muestra que para vigas muy esbeltas λ → ∞ el valor de tiende a cero. Este comportamiento nos dice que el elemento finito viga de Timoshenko de dos nodos es incapaz de reproducir en el límite la solución de la teoría clásica de vigas. Existe un fenómeno de sobre rigidización numérica que llega a “ bloquear ” la solución
haciéndola en el límite infinitamente rígida. Este elemento sólo “funciona” para vigas de relación canto/longitud elevadas y aún así su precisión no es buena, lo que lo hace inutilizable para la mayoría de los casos. Para sortear este problema, buscaremos disminuir la influencia de la cortante subintegrando los términos de K (ce ) utilizando un número de puntos de integración inferior al necesario para su cálculo exacto. Este procedimiento parece que intuitivamente realizamos una subevaluación de los términos de rigidez que tendrán como consecuencia que la flexibilidad de la viga aumente. Integraremos ahora K (ce ) con un punto de integración
(e )
Kc
⎛ GA∗ ⎞ ⎟⎟ = ⎜⎜ L ⎝ ⎠
(e )
⎡ ⎢1 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
(e )
L
−1
2
( L )
(e ) 2
4
−
(e )
L
2 1
Simétrico
⎤ ⎥ 2 ⎥ ( L(e) )2 ⎥ 4 ⎥ e L( ) ⎥ ⎥ − 2 ⎥ ( L(e) )2 ⎥ ⎥ 4 ⎦ (e )
L
(1.59)
Utilizando esta matriz de rigidez de cortante encontraremos la matriz de rigidez y flexibilidad de la viga en voladizo estudiado en este ejemplo
⎡ GA∗ ⎢ L K=⎢ ∗ ⎢− GA ⎢⎣ L
⎤ − ⎥ 2 ⎥ ⎛ GA∗ EI ⎞⎥ ⎜⎜ ⎟⎟ L + 4 L ⎝ ⎠⎥⎦ GA
∗
3 2 ⎡⎛ L L ⎞ L ⎤ ⎟⎟ ⎢⎜⎜ ∗ + ⎥ 4 2 EI EI GA ⎝ ⎠ ⎢ ⎥ y F= 2 ⎢ L L ⎥ ⎢⎣ ⎥ 2 EI EI ⎦
(1.60)
Observamos que F coincide con la expresión (1.56) a excepción del coeficiente F 11 . Resolviendo obtenemos el valor de la flecha w2
3 ⎛ L L ⎞ ⎟⎟ P = F 11 P = ⎜⎜ ∗ + 4 GA EI ⎝ ⎠
(1.61)
La relación entre este valor y el exacto para vigas esbeltas nos da: ϕ =
w2
(w2 ) f exacta
3λ 2 + 3 = 4λ 2
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
(1.62)
23
Pág.24
Elemento Finito Viga
Valores de Número de elementos ϕ =
w
L h
→∞
1
2
3
4
16
0,75
0,938
0,984
0,996
0,999
w f
Tabla 4.2.1 Convergencia con el número de elementos de la relación ϕ entre la flechas en el extremo de una viga empotrada obtenidas con y sin inclusión del efecto del esfuerzo cortante
La tabla anterior nos indica por consiguiente que la integración reducida nos proporciona un elemento válido para vigas de pequeño y gran canto. Una vez calculados los desplazamientos nodales, los esfuerzos se obtienen en el punto de Gausscentral que además es en este caso el punto óptimo.|
Grupo de Métodos Computacionales para Ingeniería GMC-PUCP
24