Análisis de Estructuras con Elementos Finitos
Análisis de Estructuras con Elementos Finitos
1.1
Cálculo de Estructuras por el Método de Elementos Finitos
1.2
SISTEMAS DISCRETOS Y CONTINUOS
(e) (e)
Finalmente, estableciendo el equilibrio de las fuerzas axiles R1 y R2 actuantes
en los extremos de la barra, se tiene (ver Figura 1.2)
(e) (e)
(e) (e) u − u1 (e) (e)
R2 = −R1 = N = (EA)(e) 2 = k(e) (u2 − u1 ) (1.4)
l(e)
(e)
donde k(e) = EA l . El ı́ndice e indica que los valores se refieren a una barra
particular. La ec.(1.4) puede escribirse en forma matricial como†
R(e)
1 −1 u1
(e)
q(e) = 1 = k(e) = K(e) a(e) (1.5)
(e)
R2 −1 1 u(e)
2
1.3
Cálculo de Estructuras por el Método de Elementos Finitos
Ka = f (1.8b)
1.4
SISTEMAS DISCRETOS Y CONTINUOS
Todos los resultados deben presentarse con la mayor claridad, y de forma gráfica
si es posible para facilitar la toma de decisiones sobre el diseño. Esta presentación
constituye la etapa de postproceso que, al igual que la de preproceso, debe estar
preparada para poder adaptarse a todas las posibles opciones de cada tipo de
problema.
1.5
Cálculo de Estructuras por el Método de Elementos Finitos
– Solución
De acuerdo con la ec.(1.5) la ecuación de equilibrio de cada barra es la siguiente:
(1)
(1)
R1 (1) 1 −1 u1
Barra 1 (1) = k (1)
R2 −1 1 u2
(2)
(2)
R1 1 −1 u1
Barra 2 (2) = k (2) (2)
R2 −1 1 u2
(3)
(3)
R1 (3) 1 −1 u1
Barra 3 (3) = k (3)
R2 −1 1 u2
EA 2EA
con k (1) = k (2) = l y k (3) = l .
1.6
SISTEMAS DISCRETOS Y CONTINUOS
(e)
Sustituyendo los valores de Ri obtenidos de las ecuaciones de equilibrio de cada
barra se llega a las ecuaciones siguientes:
(1) (1)
nudo 1 : k (1)(u1 − u2 ) = −R1
(2) (2)
nudo 2 : k (2)(u1 − u2 ) = −R2
(1) (1) (2) (2) (3) (3)
nudo 3 : k (1)(−u1 + u2 ) + k (2)(−u1 + u2 ) + k (3)(u1 − u2 ) = 0
(3) (1)
nudo 4 : k (3)(−u1 + u2 ) = P
1 2 3 4
1 k (1) 0 −k (1) 0
u1
−R1
(2) (2)
2 0 k −k 0 u2 −R2
−k (1) −k (2) (k (1) + k (2) + k (3) ) −k (3)
=
3
u3
0
4 0 0 −k (3) k (3) u4 P
Sustituyendo los valores de las rigideces de cada barra k (e) e imponiendo las
condiciones de contorno u1 = u2 = 0 se encuentra, resolviendo el sistema anterior
Pl Pl P
u3 = ; u4 = ; R1 = R2 =
2EA EA 2
y los esfuerzos axiles en cada barra
EA P
Barra 1 : N (1) = (u3 − u1 ) =
l 2
EA P
Barra 2 : N (2) = (u3 − u2 ) =
l 2
2EA
Barra 3 : N (3) = (u4 − u3 ) = P
l
1.7
Cálculo de Estructuras por el Método de Elementos Finitos
1.8
SISTEMAS DISCRETOS Y CONTINUOS
(e) (e)
donde δu1 y δu2 son, respectivamente, los desplazamientos virtuales de los
extremos 1 y 2 de la barra de volumen V (e) , y δε la correspondiente deformación
(e) (e)
virtual que puede calcularse en función de δu1 y δu2 por (1.1) como
(e) (e)
δu2 − δu1
δε = (1.14)
l(e)
que son las relaciones de equilibrio buscadas entre las fuerzas y desplazamientos
de los extremos de la barra. Como puede apreciarse, dichas ecuaciones, escritas
en forma matricial, coinciden con las (1.4) obtenidas de manera directa.
El PTV se utilizará constantemente a lo largo del libro para obtener las
ecuaciones matriciales de equilibrio de los diferentes elementos finitos para cada
tipologı́a estructural.
1.9
Cálculo de Estructuras por el Método de Elementos Finitos
Figura 1.6 Fuerzas y desplazamientos en los nudos de una barra de una estructura
articulada plana.
Si consideramos una barra 1-2 inclinada con respecto al eje global x, se deduce
para el nudo 1 que (Figura 1.6)
(e) (e) (e) (e)
Rx1 = R1 cos α ; Ry1 = R1 sen α
(e) (e) (e)
u1 = u1 cos α + v1 sen α (1.18)
donde las primas indican componentes en la dirección del eje local de la barra x .
En forma matricial
R(e)
cos α
R1 = [L(e) ]T R1
(e) x1
q1 = =
(e) senα
Ry1
(e) u1 (e) (e)
u1 = [cos α, sen α] = L(e) u1 (1.19)
v1
(e) (e)
donde u1 y q1 contienen los dos desplazamientos y fuerzas en el nudo 1 según
las direcciones cartesianas globales x e y, respectivamente, y L(e) = [cos α, sen α].
1.10
SISTEMAS DISCRETOS Y CONTINUOS
(e)
q2 = [L(e) ]T R2
(e) (e)
y u2 = L(e) u2 (1.20)
con
(e) (e) (e) (e) EA (e)
R1 = −R2 = k(e) u1 − u2 con k(e) = (1.21)
l
o, en forma matricial
(e) (e) (e) (e)
q1 K11 K12
1
u
=
(1.23)
(e)
(e) (e) (e)
q2 K21 K22 u2
en la que
1.11
Cálculo de Estructuras por el Método de Elementos Finitos
(e) (e)
K11 K12
K(e) = (e) (e)
K21 K22
i m
ui Rxi
i v R
(e) (e)
K11 K12 i yi
=
um Rxm
(e) (e)
m K21 K22 vm Rym
1.12
SISTEMAS DISCRETOS Y CONTINUOS
1 2 2 3
(1) (1) (2) (2)
K11 K12 1 K11 K12 2
K(1) = ; K(2) =
(1) (1) (2) (2)
K21 K22 2 K21 K22 3
1 2 3
(1) (1)
1 K11 K12 0 a1
f1
(e)
− − −
− − −
(1) (1) (2)
(2)
Ka = 2 K K22 + K11 K12 a2 = (e) =f
21 f2
(2)
− − −
− − −
3 0 K21
(2)
K22 a3
(e)
f 3
(e) (e)
ai = [ui, vi ]T , fi = [Rxi , Ryi ]T , Kij como en ec.(1.24)
(e) (e)
orientadas como se muestra en la Figura 1.9, y mi y θi el momento y el giro
del nudo (tomados positivos en sentido antihorario).
La deformación axial de la barra es idéntica al caso de la barra articulada
y viene definida por la ec.(1.1). Las restantes relaciones entre los esfuerzos en
los extremos y los correspondientes desplazamientos se obtienen de las ecuaciones
elásticas de la barra bajo la hipótesis de pequeños desplazamientos, que son [T7]
(e) (e)
(e) (e) (e) 3(v1 − v2 )
m1 = 2k(e) 2θ1 + θ2 +
l(e)
(e) (e)
(e) (e) (e) 3(v1 − v2 )
m2 = 2k(e) 2θ2 + θ1 + (1.26)
l(e)
1.13
Cálculo de Estructuras por el Método de Elementos Finitos
(e) (e)
(e) (e) (m1 + m2 )
Ry = −Ry = =
1 2 l(e)
(e)
12EI (e) (e) 6EI (e) (e) (e)
= 3 (v1 − v2 ) + 2 (θ1 + θ2 ) (1.27)
l l
1.14
SISTEMAS DISCRETOS Y CONTINUOS
donde
(e)
y Li es la matriz de transformación de fuerzas y movimientos globales a locales
(e) (e)
del nudo i. Debido a que la barra es recta, Li = Lj = L(e) , con (ver Figura
1.9)
cos α sen α 0
L(e) = −sen α cos α 0 (1.32)
0 0 1
De las ecs.(1.28) y (1.30) se deduce
!
[L(e) ]T 0 T
q(e) = q(e) = T(e) K(e) u(e) =
0 [L(e) ]T
T
= T(e) K(e) T(e) u(e) = K(e) u(e) (1.33)
donde !
L(e) 0
T(e) = (1.34)
0 L(e)
y
T
K(e) = T(e) K(e) T(e) (1.35)
1.15
Cálculo de Estructuras por el Método de Elementos Finitos
u2 = u2 (1.39)
1.16
SISTEMAS DISCRETOS Y CONTINUOS
De esta manera, la segunda ecuación, al ser 1015 k22 mucho mayor que el resto
de los coeficientes, equivale a
1.17
Cálculo de Estructuras por el Método de Elementos Finitos
acciones exteriores más diversas. La gran analogı́a existente entre los conceptos
del análisis matricial de estructuras de barras y los del método de los elementos
finitos facilitan en gran manera el estudio de éste a los técnicos con dominio de las
ideas sobre cálculo matricial de estructuras tratadas en apartados anteriores.
Es importante destacar desde un principio las analogı́as entre las etapas básicas
del análisis matricial de estructuras de barras y el de una estructura cualquiera
por el método de los elementos finitos. Dichas analogı́as se evidencian claramente
considerando un ejemplo, como el análisis del puente de la Figura 1.11 por
elementos finitos. Sin entrar en excesivos detalles, las etapas básicas de dicho
análisis serı́an las siguientes:
Etapa 1 : A partir de la realidad fı́sica del puente, sus apoyos y tipos de
cargas que sobre él actúen, es necesario primeramente seleccionar un modelo
matemático apropiado para describir el comportamiento de la estructura. Por
ejemplo, podrı́a utilizarse la teorı́a de láminas planas, láminas curvas, o la de
elasticidad tridimensional. También hay que definir con detalle las propiedades
mecánicas de los materiales del puente y el carácter de la deformación del mismo
(pequeños o grandes movimientos, análisis estático o dinámico, etc.). En este
curso estudiaremos únicamente problemas de equilibrio estático de estructuras
con pequeños desplazamientos y comportamiento elástico lineal de los materiales.
Asimismo, para el planteamiento de las ecuaciones de equilibrio haremos uso
siempre del Principio de los Trabajos Virtuales (PTV).
1.18
SISTEMAS DISCRETOS Y CONTINUOS
1.19
Cálculo de Estructuras por el Método de Elementos Finitos
1.20
SISTEMAS DISCRETOS Y CONTINUOS
Figura 1.12 Organigrama general del análisis de una estructura por el método de
los elementos finitos.
1.21
Cálculo de Estructuras por el Método de Elementos Finitos
1.22
CAPÍTULO 2
ELEMENTOS FINITOS DE
BARRA. CONCEPTOS BÁSICOS
2.1 INTRODUCCIÓN
2.1
Cálculo de Estructuras por el Método de Elementos Finitos
2.2
ELEMENTOS FINITOS DE BARRA
2.3
Cálculo de Estructuras por el Método de Elementos Finitos
(1) (1)
Lógicamente u(x) tiene que tomar en los nodos 1 y 2 los valores u1 y u2 ,
respectivamente. Es decir
(1) (1) (1) (1)
u(x1 ) = u1 y u(x2 ) = u2 (2.8)
(1) (1)
siendo x1 y x2 las coordenadas de los nodos 1 y 2. El ı́ndice 1 indica que los
valores se refieren al elemento número 1.
(1) (1)
u1 = ao + a1 x1
(2.9)
(1) (1)
u2 = ao + a1 x2
2.4
ELEMENTOS FINITOS DE BARRA
2.5
Cálculo de Estructuras por el Método de Elementos Finitos
Obsérvese que por ser las funciones de forma lineales la deformación es constante
sobre el elemento.
Las fuerzas entre elementos se transmiten únicamente a través de los nodos.
Dichas fuerzas nodales, que denominaremos “de equilibrio”, pueden calcularse para
cada elemento haciendo uso del PTV, que se escribe para el elemento considerado
como
x(1) x(1)
2 2 (1) (1) (1) (1)
(1) δε EAε dx = (1) δub dx + δu1 X1 + δu2 X2 (2.15)
x1 x1
x(1)
2 (1) (1) (1) (1) (1) (1) (1) (1)
− (1) N1 δu1 + N2 δu2 b dx = δu1 X1 + δu2 X2 (2.18)
x1
y, agrupando términos
(2.19)
2.6
ELEMENTOS FINITOS DE BARRA
(1) (1)
Del sistema de ecuaciones anterior se deducen los valores de X1 y X2 . En
forma matricial
dN (1) dN1
(1) dN (1) (1)
x(1) dN
u
1 1 2 (1)
2 (EA) dx (EA) dx
dx dx dx
(1) −
1
(1) (1)
dN1 dN (1)
(1) (1)
x1 dN2 2 dN 2 u2
dx (EA) dx dx (EA) dx (2.21)
x(1) (1) X (1)
2 N1
− (1) b dx = 1
N2 (1)
x1 (1)
X2
con
x(1)
(1) 2 dNi (1) dNj (1)
Kij = (1) (EA) dx (2.22)
x1 dx dx
x(1)
(1) 2 (1)
fi = (1) Ni b dx i, j = 1, 2
x1
(1) (1) T (1) (2) T
a(1) = [u1 , u2 ] ; q(1) = [X1 , X2 ]
donde K(1) , a(1) , f (1) y q(1) son la matriz de rigidez, el vector de desplazamientos
nodales, el vector de fuerzas nodales equivalentes y el vector de fuerzas nodales de
equilibrio del elemento, respectivamente.
2.7
Cálculo de Estructuras por el Método de Elementos Finitos
expresiones que coinciden con las obtenidas para la barra bajo cargas axiles en
el Capı́tulo 1. Dicha coincidencia no es fortuita, y podı́a haberse anticipado,
ya que en ambos casos se parte de la misma hipótesis de distribución lineal de
desplazamientos, lo que evidentemente conduce a idénticas expresiones para la
matriz de rigidez y el vector de fuerzas en los extremos de la barra.
Las ecuaciones que expresan el equilibrio global de la estructura se pueden
obtener por un proceso idéntico al explicado para las estructuras de barras en
el Capı́tulo 1. Ası́, en cada nodo se tiene que satisfacer la ecuación básica de
equilibrio de fuerzas
(e)
Xi = Xjext (2.24)
e
donde el sumatorio se extiende sobre todos los elementos que concurren en el nodo
(e)
en cuestión, Xi es la fuerza de equilibrio que aporta cada elemento y Xjext la
fuerza puntual exterior sobre el nodo de número global j.
Para la malla de un solo elemento que se considera, la ec.(2.24) se escribe,
teniendo en cuenta la Figura 2.2, como
(1)
nodo 1 : X1 = R1
(1)
nodo 2 : X2 = P
l bl
u2 = (P + ) ; R1 = −(P + bl) (2.26)
EA 2
2.8
ELEMENTOS FINITOS DE BARRA
bl
N (1) = (EA)(1) ε(1) = P + (2.27)
2
2.9
Cálculo de Estructuras por el Método de Elementos Finitos
Elemento 1 Elemento 2
(1) (1) (1) (1) (2) (2) (2) (2)
u(x) = N1 (x)u1 + N2 (x)u2 u(x) = N1 (x)u1 + N2 (x)u2 (2.29)
2.10
ELEMENTOS FINITOS DE BARRA
Elemento 1 Elemento 2
Las funciones de forma y sus derivadas son ahora
(1) (1) (2) (2)
(1) x2 − x dN1 1 (2) x2 − x dN1 1
N1 = ; = − (1) N1 = ; = − (2)
l(1) dx l l(2) dx l (2.30)
(1) (1) (2) (2)
(1) x − x1 dN2 1 (2) x − x1 dN2 1
N2 = ; = (1) N2 = ; = (2)
l(1) dx l l(2) dx l
La deformación axial en un punto cualquiera de cada elemento es
(1) (1) (2) (2)
dN1 (1) dN2 (1) dN1 (2) dN2 (2)
ε = du = u + u ε = du = u + u (2.31)
dx dx 1 dx 2 dx dx 1 dx 2
2.11
Cálculo de Estructuras por el Método de Elementos Finitos
# $(1) # EA $(1)
EA
# $−(1) l # $(2) 0 bl
l
# EA $(1) # $(2) u1 4 + R1
− EA
+ EA − EA u2 = bl
(2.35)
l l l l 2
# EA $(2) # EA $(2) u3 bl
+P
0 − l l
4
Ka = f (2.36)
2.12
ELEMENTOS FINITOS DE BARRA
! "
EA (e) 1 −1
= (2.37)
l −1 1
y ensamblando seguidamente las matrices individuales de todos los elementos
siguiendo precisamente las mismas reglas del Capı́tulo 1 para las estructuras de
barras.
El mismo proceso es aplicable al vector de fuerzas nodales equivalentes.
Por tanto, si sobre los elementos actuan fuerzas uniformemente repartidas, el
ensamblaje del vector f puede efectuarse a partir del vector de fuerzas nodales
equivalentes de los diferentes elementos dado por
% &(e)
f (e) x(e) (e)
2 N1 bl 1
f (e) = 1 = (e) b(e) dx = (2.38)
(e)
N2
x1 (e) 2 1
f2
(1)
Sustituyendo EA = EA (2) = 2EA y resolviendo el sistema con u = 0
l l l 1
se encuentra
% &
l 3bl
u1 = 0 ; u2 = P+
2EA 4
(2.39)
l
u3 = (2P + bl) ; R1 = −(P + bl)
2EA
La deformación y el esfuerzo axil en cada elemento se obtienen por
Elemento 1 Elemento 2
2.13
Cálculo de Estructuras por el Método de Elementos Finitos
2.14
ELEMENTOS FINITOS DE BARRA
u = N1 u1 + N2 u2 (2.45)
2.15
Cálculo de Estructuras por el Método de Elementos Finitos
donde
u1
N = [N1, N2 ] ; a(e) = (2.47)
u2
donde
dN1 dN2
B = [ , ] (2.49)
dx dx
donde
D = [EA] (2.51)
σ = D · B · a(e) (2.52)
t×1 t × t [t × (n × d)] [(n × d) × 1]
2.16
ELEMENTOS FINITOS DE BARRA
donde
K(e) = (e)
BT D B dx
l (2.58)
T
f (e) = N b dx
l(e)
2.17
Cálculo de Estructuras por el Método de Elementos Finitos
D = [EA] y b = {b}
y sustituyendo en (2.58)
1 EA (e)
(e)
1
− l(e) 1 1 −1
K = 1 (EA) − (e) , (e) dx =
l(e) l(e) l l l −1 1
(2.60)
(e)
x2 − x b (bl)(e) 1
f (e) = dx =
l(e)
x − x1 l 2 1
T
.. T
B1T
B1 DB1 . B1 DB2
K(e) = D [B1, B2 ] dx = ......... ... . . . . . . . . . dx
l(e)
BT2 l(e) ..
BT2 DB1 . BT2 DB2 (2.61)
NT1 NT1 b
f (e) = b dx = dx
l(e)
NT2 l(e)
NT2 b
(e)
De (2.61) puede definirse la matriz Kij que relaciona los nodos i y j del
elemento como
(e)
Kij = BTi D Bj dx ; i, j = 1, 2 (2.62)
l(e)
d.d (d × t) (t × t) (t × d)
2.18
ELEMENTOS FINITOS DE BARRA
(e)
fi = NTi b dx i = 1, 2 (2.63)
l(e)
(d × 1) (d × d) (d × 1)
Ka = f
2.19
Cálculo de Estructuras por el Método de Elementos Finitos
a = K−1f (2.67)
2.20
CAPÍTULO 3
ELEMENTOS DE BARRA
MÁS AVANZADOS
3.1 INTRODUCCIÓN
3.1
Cálculo de Estructuras por el Método de Elementos Finitos
u(x) = αo + α1 x + α2 x2 + · · · (3.1)
u(x1) = u1 = αo + α1 x1
(3.3)
u(x2) = u2 = αo + α1 x2
donde u1 y u2 , son los valores del desplazamiento axial en los nodos. Despejando
αo y α1 y sustituyendo en (3.1), se obtiene
3.2
ELEMENTOS DE BARRA MÁS AVANZADOS
donde
(x2 − x) (x − x1 )
N1 (x) = ; N2 (x) = (3.5)
l(e) l(e)
x − x2 x −x x − x1 x − x1
N1 = = 2(e) , N2 = = (3.7)
x1 − x2 l x2 − x1 l(e)
3.3
Cálculo de Estructuras por el Método de Elementos Finitos
(ξ − ξ2 )(ξ − ξ3 ) 1
N1 = = ξ(ξ − 1)
(ξ1 − ξ2 )(ξ1 − ξ3 ) 2
(ξ − ξ1 )(ξ − ξ3 )
N2 = = (1 + ξ)(1 − ξ) (3.11)
(ξ2 − ξ1 )(ξ2 − ξ3 )
(ξ − ξ1 )(ξ − ξ2 ) 1
N3 = = ξ(1 + ξ)
(ξ3 − ξ1 )(ξ3 − ξ2 ) 2
(ξ − ξ2 )(ξ − ξ3 )(ξ − ξ4 ) 9 1 1
N1 = = − (ξ + )(ξ − )(ξ − 1)
(ξ1 − ξ2 )(ξ1 − ξ3 )(ξ1 − ξ4 ) 16 3 3
(ξ − ξ1 )(ξ − ξ3 )(ξ − ξ4 ) 27 1
N2 = = (ξ + 1)(ξ − )(ξ − 1)
(ξ2 − ξ1 )(ξ2 − ξ3 )(ξ2 − ξ4 ) 16 3
(3.12)
(ξ − ξ1 )(ξ − ξ2 )(ξ − ξ4 ) 27 1
N3 = = − (ξ + 1)(ξ + )(ξ − 1)
(ξ3 − ξ1 )(ξ3 − ξ2 )(ξ3 − ξ4 ) 16 3
(ξ − ξ1 )(ξ − ξ2 )(ξ − ξ3 ) 9 1 1
N4 = = (ξ + 1)(ξ + )(ξ − )
(ξ4 − ξ1 )(ξ4 − ξ2 )(ξ4 − ξ3 ) 16 3 3
3.4
ELEMENTOS DE BARRA MÁS AVANZADOS
3.3.1 Introducción
Una vez estudiada la obtención general de las funciones de forma de los
elementos unidimensionales de clase Co más usuales, es el momento oportuno
de introducir dos importantes conceptos sin los cuales es prácticamente imposible
que el método de los elementos finitos se hubiese desarrollado hasta los niveles en
que hoy se encuentra.
El primer concepto es el de formulación paramétrica. La idea es interpolar
la geometrı́a del elemento a partir de las coordenadas de una serie de puntos
3.5
Cálculo de Estructuras por el Método de Elementos Finitos
3.6
ELEMENTOS DE BARRA MÁS AVANZADOS
punto del mismo interpolando los valores de las coordenadas conocidas. Dicha
interpolación puede escribirse en la forma
3.7
Cálculo de Estructuras por el Método de Elementos Finitos
con lo que
l(e) dξ 2
dx = dξ y = (e) (3.20)
2 dx l
dN1 2 dN1 1
= (e) = − (e)
dx l dξ l (3.21)
dN2 2 dN2 1
= (e) = (e)
dx l dξ l
y, por consiguiente
1 1
B = − (e) , (e) (3.22)
l l
Resultado que, por otra parte, ya conocı́amos. Hay que resaltar que en este caso
particular pueden obtenerse las ecs.(3.20) de una manera más sencilla a partir de
(3.8). No obstante, hemos preferido seguir aquı́ un procedimiento más sistemático
que facilitará la comprensión del desarrollo de elementos isoparamétricos más
complejos.
La matriz de rigidez y el vector de fuerzas nodales equivalentes de la ec.(2.58)
se pueden expresar ahora en el sistema normalizado haciendo uso de (3.20) como
+1
l(e)
K(e) = BT (EA) B dξ
−1 2
+1 (3.23)
l(e)
f (e) = NT b dξ
−1 2
3.8
ELEMENTOS DE BARRA MÁS AVANZADOS
3 dN
u1
du dNi 1dξ dN2 dξ dN3 dξ
ε = = ui = , , ]
u 2
= B a(e)
u3
dx i=1 dx dξ dx dξ dx dξ dx
(3.26)
Por otra parte, de (3.11) se deduce
dN1 1 dN2 dN3 1
= ξ− ; = −2ξ ; =ξ+ (3.27)
dξ 2 dξ dξ 2
dξ 2
= (3.30)
dx l + 2ξ(x1 + x3 − 2x2)
(e)
3.9
Cálculo de Estructuras por el Método de Elementos Finitos
parte usual) de que el nodo intermedio esté situado en el centro del elemento, se
tiene
dξ 2
= (3.31)
dx l (e)
y, por consiguiente
dx l(e) l(e)
= y dx = dξ (3.32)
dξ 2 2
2 1 1
B = (ξ − ), −2ξ, (ξ + ) (3.33)
l(e) 2 2
+1 − 12 )
(ξ
2 2 1 1 l(e)
K(e) = −2ξ (EA) [(ξ − ), −2ξ, (ξ + )] dξ (3.35)
−1 l(e)
1
l (e) 2 2 2
(ξ + 2 )
3.10
ELEMENTOS DE BARRA MÁS AVANZADOS
3.11
Cálculo de Estructuras por el Método de Elementos Finitos
n ±ξi Wi
1 0.0 2.0
2 0.5773502692 1.0
0.774596697 0.5555555556
3 0.0 0.8888888889
0.8611363116 0.3478548451
4 0.3399810436 0.6521451549
0.9061798459 0.2369268851
5 0.5384693101 0.4786286705
0.0 0.5688888889
0.9324695142 0.1713244924
6 0.6612093865 0.3607615730
0.2386191861 0.4679139346
0.9491079123 0.1294849662
0.7415311856 0.2797053915
7 0.4058451514 0.3818300505
0.0 0.4179591837
0.9602898565 0.1012285363
0.7966664774 0.2223810345
8 0.5255324099 0.3137066459
0.1834346425 0.3626837834
3.12
ELEMENTOS DE BARRA MÁS AVANZADOS
Figura 3.3 Puntos óptimos para cálculo de tensiones en algunos elementos uni y
bidimensionales.
3.13
Cálculo de Estructuras por el Método de Elementos Finitos
3.14
ELEMENTOS DE BARRA MÁS AVANZADOS
Figura 3.4 Ejemplo de aproximación de una solución cúbica con diferentes tipos
de elementos finitos. Para mayor sencillez se ha supuesto que en todos
los casos la aproximación es exacta en los nodos.
u = N1 u1 + N2 u2 + . . . + Nn un =
u 1
n u2 (3.41)
= Ni ui = [N1 , N2 , . . . , Nn ] .. = N a(e)
.
i=1
un
3.15
Cálculo de Estructuras por el Método de Elementos Finitos
Por otra parte, las derivadas cartesianas de las funciones de forma se obtienen
por
dNi dNi dξ
= (3.44)
dx dξ dx
De (3.42) se deduce
dx n dN
i
= xi = J (e) (3.45)
dξ i=1 dξ
Por consiguiente
dξ 1
dx = J (e) dξ ; = (e) (3.46)
dx J
y
dNi 1 dNi
= (e) (3.47)
dx J dξ
3.16
ELEMENTOS DE BARRA MÁS AVANZADOS
con D = [EA]
3.17
Cálculo de Estructuras por el Método de Elementos Finitos
3.18
ELEMENTOS DE BARRA MÁS AVANZADOS
Definición de:
• Tipo de elemento
SUBRUTINA DATOS
• Topologı́a y coordenadas de la malla
Entrada • Propiedades del material
• Condiciones de contorno
de datos geométricos
y del material. • Coordenadas y pesos de la
cuadratura de Gauss-Legendre
Cálculo en cada punto de Gauss de:
• Propiedades del material(EA)
• Derivadas ∂N i
SUBRUTINA RIGIDEZ
∂ξ
∂Ni
Cálculo de la matriz • J (e) = ∂ξ i
x
i
de rigidez de cada
• Matriz B
elemento.
Obtenciónde:
(e)
[J (e) BT (EA)B]p Wp
• K =
p
Cálculo en cada punto de Gauss de:
• Funciones de forma Ni
• J (e) = ∂Ni x
SUBRUTINA CARGAS
Cálculo del vector de ∂ξ i
i
fuerzas nodales
Obtención
de:(e) T
• f (e) = [J N b]p Wp
de cada elemento. p
SUBRUTINA SOLUCION
Eliminación Gaussiana [R2]
Ensamblaje y Método Frontal [H5]
solución del sistema. Método del perfil [Z6], etc.
Ka=f
SUBRUTINA TENSION
Cálculo de deformaciones ε = Ba
σ = DBa
y tensiones en cada
elemento.
STOP
3.19
CAPÍTULO 4
SÓLIDOS BIDIMENSIONALES
4.1 INTRODUCCIÓN
4.1
Cálculo de Estructuras por el Método de Elementos Finitos
4.2
SÓLIDOS BIDIMENSIONALES
donde u(x, y) y v(x, y) son los desplazamientos del punto en direcciones de los ejes
x e y, respectivamente.
4.3
Cálculo de Estructuras por el Método de Elementos Finitos
σ=Dε (4.5)
E E(1 − ν)
d11 = d22 = d11 = d22 =
1 − ν2 (1 + ν)(1 − 2ν)
ν
d12 = d21 = νd11 d12 = d21 = d11 (4.7)
1−ν
E E
d33 = =G d33 = =G
2(1 + ν) 2(1 + ν)
4.4
SÓLIDOS BIDIMENSIONALES
σ = D εe = D (εε − ε0 ) (4.8)
Para el caso usual de deformación inicial isótropa por efectos térmicos el vector
0
ε tiene la expresión siguiente:
donde
4.5
Cálculo de Estructuras por el Método de Elementos Finitos
(4.12)
+ (δutx + δvty )t ds + (δui Ui + δvi Vi )
l i
El segundo miembro representa el trabajo de las fuerzas repartidas por unidad
de volumen bx , by ; de las fuerzas repartidas sobre el contorno tx , ty ; y de las fuerzas
puntuales Ui , Vi sobre los desplazamientos virtuales δu, δv. El primer miembro,
por otro lado, representa el trabajo que las tensiones σx , σy , τxy realizan sobre
las deformaciones virtuales δεx , δεy y δγxy . A y l son el área y el contorno de
la sección transversal del sólido y t su espesor. En problemas de tensión plana t
coincide con el espesor real, mientras que en problemas de deformación plana es
usual asignar a t un valor unidad.
La ec.(4.12) se puede reescribir en forma matricial como
δεεT σ
ε σ tdA = δuT b tdA + δuT t tds + δuTi qi (4.13)
A A l i
donde
T T T
δεε = δεx , δεy , δγxy ; δu = δu, δv ; b = bx , by
T T T
t = tx , ty ; δui = δui , δvi ; qi = Ui , Vi
(4.14)
4.6
SÓLIDOS BIDIMENSIONALES
4.7
Cálculo de Estructuras por el Método de Elementos Finitos
u = N1 u1 + N2 u2 + N3 u3
(4.15)
v = N1 v1 + N2 v2 + N3 v3
4.8
SÓLIDOS BIDIMENSIONALES
donde
u
u= (4.18)
v
son la matriz de funciones de forma del elemento y del nodo i del elemento,
respectivamente, y
(e)
a
1
ui
(e)
a(e) = a
(e)
2
con ai = (4.20)
(e)
vi
a3
u = α1 + α2 x + α3 y
(4.21)
v = α4 + α5 x + α6 y
u1 = α1 + α2 x1 + α3 y1
u2 = α1 + α2 x2 + α3 y2 (4.22)
u3 = α1 + α2 x3 + α3 y3
1
u= (a 1 + b 1 x + c1 y)u 1 + (a 2 + b 2 x + c2 y)u 2 + (a 3 + b 3 x + c3 y)u 3 (4.23)
2A(e)
4.9
Cálculo de Estructuras por el Método de Elementos Finitos
ai = xj yk − xk yj , bi = yj − yk , ci = xk − xj ; i, j, k = 1, 2, 3 (4.24)
Comparando (4.23) con (4.15) se deduce que las funciones de forma del elemento
son
1
Ni = (ai + bi x + ci y) , i = 1, 2, 3 (4.25)
2A(e)
La representación gráfica de dichas funciones se muestra en la Figura 4.4. Puede
comprobarse, como ejercicio, que las funciones de forma toman el valor unidad en
un nodo y cero en los otros dos.
4.10
SÓLIDOS BIDIMENSIONALES
y en forma matricial
u1
.. ..
∂u
∂N1
∂x 0 . ∂N2 0 . ∂N3 0
v1
∂x
∂x ∂x
∂v ∂N1 .. ∂N2 .. ∂N3 u2
ε= ∂y == 0 . 0 . 0 (4.27)
∂y ∂y ∂y v2
∂u + ∂v
∂N1 ∂N1 ..
. ∂N2 ∂N2 ..
. ∂N3 ∂N3
u3
∂y ∂x ∂y ∂x ∂y ∂x ∂y ∂x v3
o
ε = Ba(e) (4.28)
donde
B = [B1 , B2 , B3 ] (4.29)
y, por consiguiente
b 0
1 i
Bi = 0 ci (4.32)
2A(e) c bi
i
σ = D ε = D B a(e) (4.33)
4.11
Cálculo de Estructuras por el Método de Elementos Finitos
4.12
SÓLIDOS BIDIMENSIONALES
3
3
(e)
δεεT σ tdA = (e)
δuT b tdA + (e)
δuT t tds + δuiUi + δviVi (4.35)
A A l i=1 i=1
donde δui y δvi son los desplazamientos virtuales de los nodos del elemento y Ui
y Vi las fuerzas nodales de equilibrio que corresponden a dichos desplazamientos.
El trabajo virtual de dichas fuerzas puede despejarse de la ecuación anterior como
T
(e)
δεεT σ t dA − (e)
δuT b t dA − (e)
δuT t t ds = [δa(e) ] q(e) (4.36)
A A l
T T
[δa(e)] (e)
BT σ tdA− (e)
NT b tdA− (e)
NT t tdS = [δa(e)] q(e) (4.39)
A A l
0 0
BT (DBa(e) − Dεε + σ ) tdA − NT b tdA − NT t tds = q(e)
A(e) A(e) l(e)
(4.41)
4.13
Cálculo de Estructuras por el Método de Elementos Finitos
y operando
(e)
BT D B t dA a(e) − BT Dεε0 tdA +
A
A(e) (4.42)
+ BT σ 0
σ tdA − NT b tdA − NT t tdS = q(e)
A(e) A(e) l (e)
o
K(e) a(e) − f (e) = q(e) (4.43)
donde
K(e) = BT D B tdA (4.44)
A(e)
(e)
fσ =− BT σ 0 tdA (4.47)
A(e)
(e)
fb = NT b tdA (4.48)
A(e)
(e)
ft = NT t tds (4.49)
l(e)
4.14
SÓLIDOS BIDIMENSIONALES
Ka=f (4.51)
(e)
Por consiguiente, una submatriz de rigidez tı́pica, Kij , que relacione los nodos
i y j del elemento se puede calcular como
(e)
Kij = BTi DBj t dA (4.53)
A(e)
4.15
Cálculo de Estructuras por el Método de Elementos Finitos
d11 d12 0
bj 0
(e) 1 bi 0 ci 1
Kij = d21 d22 0 0 cj tdA (4.54)
A(e) 2A (e) 0 ci bi
0 0 d33 2A(e) cj bj
4.16
SÓLIDOS BIDIMENSIONALES
(e)
En el cálculo de ft hay que tener en cuenta que al referirse la integral a un
lado del elemento, la función de forma del nodo no perteneciente a dicho lado vale
cero sobre el mismo. Ası́, si el lado cargado es el 1-2 y las fuerzas tx y ty están
uniformemente repartidas sobre dicho lado, se obtiene de (4.60) que la fuerza total
sobre el lado se reparte equitativamente entre los dos nodos del mismo y el vector
(e)
ft es
tx
ty
(e)
(e) (l12 t) tx
ft = (4.61)
2
ty
0
0
(e)
donde l12 es la longitud del lado 1-2. Se deduce fácilmente que si los lados cargados
(e)
son el 1-3 y el 2-3, la expresión de ft es en cada caso
tx
0
ty
0
(e)
(e)
(e) (l13 t) 0 (e) (l23 t) tx
ft = ; ft = (4.62)
2
0 2
ty
t
t
x
x
ty ty
t(e) bi (d11 ε0x + d12 ε0y ) + ci d33 γxy
0
= (4.65)
2 ci (d21 ε0x + d22 ε0y ) + bi d33 γxy
0
4.17
Cálculo de Estructuras por el Método de Elementos Finitos
4.18
SÓLIDOS BIDIMENSIONALES
dξ 1 dη 1
= ; = (4.73)
dx a dy b
dx dy = ab dξ dη (4.74)
4.19
Cálculo de Estructuras por el Método de Elementos Finitos
Por tanto, para integrar una función f (x, y) sobre un elemento rectangular,
puede efectuarse la siguiente transformación al sistema de coordenadas naturales
+1 +1
f (x, y)dx dy = g(ξ, η)ab dξ dη (4.75)
A(e) −1 −1
4.20
SÓLIDOS BIDIMENSIONALES
η, como se indica en la Figura 4.8. Se observa que las funciones de forma no son
nunca polinomios completos y todas contienen un número de términos adicionales
que crece con el orden del elemento.
4.21
Cálculo de Estructuras por el Método de Elementos Finitos
1 1
l1i (ξ) = (1 + ξξi ) ; l1i (η) = (1 + ηηi ) (4.77)
2 2
donde ξi y ηi toman los valores de la tabla de la Figura 4.9. Por consiguiente, la
función de forma del nodo i es
1
Ni (ξ, η) = l1i (ξ)l1i (η) = (1 + ξξi )(1 + ηηi ) (4.78)
4
En la Figura 4.9 se muestra de forma gráfica la obtención de la función de forma
del nodo 1.
4.22
SÓLIDOS BIDIMENSIONALES
a) Nodos esquina
1
Ni = (ξ 2 + ξξi )(η2 + ηηi ) ; i = 1, 3, 5, 7 (4.81)
4
1 1
Ni = ηi2 (η2 − ηηi )(1 − ξ 2 ) + ξi2 (ξ 2 − ξξi )(1 − η2 ) ; i = 2, 4, 6, 8 (4.82)
2 2
c) Nodo central
4.23
Cálculo de Estructuras por el Método de Elementos Finitos
4.24
SÓLIDOS BIDIMENSIONALES
del mismo grado que la variación sobre los lados. En la Figura 4.11 se muestran
algunos de los elementos rectangulares de la familia Serendı́pita más populares,
ası́ como los términos que intervienen en sus funciones de forma. Se observa que el
elemento más sencillo de esta familia es el rectángulo de cuatro nodos ya estudiado
y que, por consiguiente, pertenece a ambas familias, Lagrangiana y Serendı́pita.
Se aprecia, asimismo, que los elementos cuadrático y cúbico, de 8 y 12 nodos,
respectivamente, no tienen nodos internos, mientras que el de 17 nodos precisa un
nodo en su interior para poder conseguir todos los términos del polinomio completo
de cuarto grado [O3].
Las caracterı́sticas de los elementos Serendı́pitos impiden que sus funciones
de forma puedan obtenerse de un modo tan sistemático como las de los elementos
Lagrangianos. De hecho, dichas funciones de forma suelen obtenerse en la práctica
combinando la observación y el ingenio. De ahı́ la denominación Serendı́pita para
esta familia de elementos como referencia a los descubrimientos ingeniosos del
prı́ncipe de Serendip, citado en los romances de Horacio Walpole en el siglo XVIII.
No obstante, para los elementos más populares de la familia Serendı́pita, que
se presentan en la Figura 4.11, la obtención de las funciones de forma es sencilla,
como inmediatamente comprobaremos para el elemento de ocho nodos.
4.25
Cálculo de Estructuras por el Método de Elementos Finitos
Etapa 2 . Se impone que la función de forma sea nula en uno de los nodos
adyacentes, restando a la ec.(4.85) la mitad del valor de la función de forma en
dicho nodo. Ası́ , para anular N1L en el nodo 2, se hace
1
N 1 (ξ, η) = N1L − N2 (4.86)
2
4.26
SÓLIDOS BIDIMENSIONALES
1 1
N1 (ξ, η) = N1L − N2 − N8 (4.87)
2 2
4.27
Cálculo de Estructuras por el Método de Elementos Finitos
Realizando las etapas anteriores para el resto de los nodos esquina puede
encontrarse la expresión general de las funciones de forma como
1
Ni (ξ, η) = (1 + ξξi )(1 + ηηi )(ξξi + ηηi − 1) i = 1, 3, 5, 7 (4.88)
4
φ = αo + α1 x + α2 y + α3 xy + α4 x2 + α5 y2 (4.89)
φ = αo + α1 x + α2 y + α3 xy + α4 x2 + α5 y2 + α6 x3 + α7 x2y + α8 xy2 + α9 y3
(4.90)
4.28
SÓLIDOS BIDIMENSIONALES
L1 + L2 + L3 = 1 (4.92)
4.29
Cálculo de Estructuras por el Método de Elementos Finitos
x = L1 x 1 + L2 x 2 + L3 x 3
(4.93)
y = L1 y1 + L2 y2 + L3 y3
donde A(e) es el área del triángulo y ai, bi y ci coinciden con los valores de
(4.25) para el elemento triangular de tres nodos. Se comprueba, por tanto, que
las coordenadas de área son precisamente las funciones de forma del elemento
triangular de tres nodos.
4.30
SÓLIDOS BIDIMENSIONALES
Con el objeto de que la ecuación (4.95) sea consistente para todos los nodos, es
necesario tener en cuenta que Ls1 0 (L1 ) = 1.
La dificultad mayor para aplicar la ec. (4.96) consiste en deducir los valores
I, J, K de cada nodo. Esto puede hacerse fácilmente teniendo en cuenta que:
a) La función de forma de cada nodo de vértice depende únicamente de una
coordenada de área, de lo que se deduce el exponente que afecta a dicha función y,
por tanto, el valor de I, J o K del nodo; b) Los nodos colocados sobre las rectas
L1 = cte tienen el mismo I, ocurriendo lo mismo con L2 y J y L3 y K, y c) Los
valores I, J, K asociados a L1 , L2 y L3 decrecen de unidad en unidad desde sus
valores máximos sobre las rectas Li = 1 que pasan sobre los nodos de vértice,
hasta el valor cero sobre la recta Li = 0 que coincide con el lado opuesto al vértice
en cuestión. (Figura 4.14).
Aclararemos todos estos conceptos con varios ejemplos.
4.31
Cálculo de Estructuras por el Método de Elementos Finitos
Nodo 1
Posición (I, J, K) : (1, 0, 0)
Coordenadas de área: (1, 0, 0)
Nodo 4
Posición (I, J, K) : (1, 1, 0)
Coordenadas de área: (1/2, 1/2, 0)
L1 − Lj1 L2 − Lj2
N4 = l12 (L1 ) l12 (L2 )l03 (L3 ) = l12 (L1 )l12 (L2 ) = =
j=1,3 L21 − Lj1 j=1,3
j
L22 − L2
j=1,2 j=1,2
L1 − L31 L2 − L32
(L1 − 0)L1 L2 − 0
= = = 4L1 L2 (4.99)
L21 − L31 L22 − L32 (1/2) − 0 1/2 − 0
4.32
SÓLIDOS BIDIMENSIONALES
4.33
Cálculo de Estructuras por el Método de Elementos Finitos
Figura 4.17 Placa traccionada por carga parabólica. Análisis con elementos
triangulares de 3 y 6 nodos y rectangulares de 4 y 8 nodos [G2], [Y1].
4.34
SÓLIDOS BIDIMENSIONALES
Figura 4.18 Viga en voladizo bajo carga parábolica en el borde. Análisis con
elementos triangulares de 3 y 6 nodos, rectangulares de 4 y 8 nodos
y el rectangular de 4 nodos y dos modos incompatibles [O3].
4.35
Cálculo de Estructuras por el Método de Elementos Finitos
n
n
x= Ni (ξ, η) xi , y= Ni (ξ, η) yi (4.102)
i=1 i=1
donde Ni (ξ, η) son precisamente las funciones de forma del elemento. Las
ecs.(4.102) relacionan las coordenadas cartesianas de un punto y las naturales ξ y η.
Dicha relación debe ser biunı́voca, para lo cual debe cumplirse que el determinante
de la matriz Jacobiano de la transformación de coordenadas xy → ξη (dicha matriz
se define más adelante) sea de signo constante en todo el elemento [S4].
Puede demostrarse que si se utilizan funciones de forma lineales dicha condición
exige que ningún ángulo interior entre dos lados del elemento sea mayor o igual
que 180◦ [S4]. Si las funciones de forma son cuadráticas es necesario además que
los nodos sobre los lados se encuentren en el tercio central de la distancia entre
los nodos esquina adyacentes [J3]. Para funciones de forma de órdenes superiores
no existen reglas prácticas y es necesario comprobar el signo del determinante del
Jacobiano. No obstante, las funciones de grado superior a dos son poco utilizables
en la práctica. En la Figura 4.19 se muestran algunos ejemplos de elementos
isoparamétricos en dos dimensiones.
Gran parte de las ideas subyacentes en la aproximación isoparamétrica son
originales de Taig [T1], [T2], quien las aplicó para obtener siempre elementos
cuadriláteros de 4 nodos. Posteriormente, Irons [I1,4] generalizó dichas ideas para
obtener otros elementos más complejos en dos y tres dimensiones.
4.36
SÓLIDOS BIDIMENSIONALES
o en forma matricial
∂Ni
∂x ∂y
∂Ni
∂Ni
∂ξ ∂ξ ∂ξ ∂x (e) ∂x
= ∂N = J (4.104)
∂Ni
∂x ∂y i
∂Ni
∂η ∂η ∂η ∂y ∂y
!
J(e)
4.37
Cálculo de Estructuras por el Método de Elementos Finitos
" "
" "
donde "J(e) " es el determinante del Jacobiano.
El determinante del Jacobiano permite también expresar el diferencial de área
en coordenadas naturales como [C9]
" "
" "
dx dy = "J(e) " dξ dη (4.106)
4.38
SÓLIDOS BIDIMENSIONALES
∂x n ∂N (α, β)
∂x n ∂N (α β)
i i 1
= xi ; = xi ; etc. (4.113)
∂α i=1 ∂α ∂β i=1 ∂β
donde todos los términos de Bi , J(e) y Gij se deducen de las ecs. (4.107) - (4.110)
sustituyendo las coordenadas ξ y η por α y β, respectivamente.
En elementos de lados curvos los términos del integrando de (4.114) son
funciones racionales en α y β. Esta dificultad se complica con la interdependencia
de los lı́mites de integración debida a la geometrı́a del elemento. No obstante, el
cálculo de las integrales puede efectuarse de manera sencilla y sistemática mediante
integración numérica.
4.39
Cálculo de Estructuras por el Método de Elementos Finitos
+1 +1 +1
nq np nq
g(ξ, η) dξdη = g(ξ, ηq )Wq dξ = g(ξp , ηq )Wp Wq
−1 −1 −1 q=1 p=1 q=1
(4.115)
donde np y nq son el número de puntos de integración seleccionados en cada una de
las direcciones ξ y η; ξp y ηq son las coordenadas naturales del punto de integración
p, q y Wp, Wq los pesos correspondientes a cada dirección en dicho punto.
Las coordenadas y los pesos para cada dirección se deducen directamente de
los dados en la Tabla 3.1 para el caso unidimensional. Recordemos que una
cuadratura de orden n en cada dirección natural integra exactamente un polinomio
de grado ≤ 2n − 1 en la correspondiente coordenada natural. En la Figura 4.20 se
muestran algunas de las cuadraturas bidimensionales más usuales sobre elementos
cuadriláteros.
1 1−L np
3
f (L1 , L2 , L3) dL2 dL3 = f (L1p , L2p , L3p ) Wp (4.116)
0 0 p=1
donde np es el número de puntos de integración; L1p , L2p , L3p y Wp son los valores
de las coordenadas de área y del peso en el punto de integración p, respectivamente.
En la Figura 4.21 se muestran las coordenadas y los pesos más utilizados en
la práctica; la “precisión”en dicha figura es el polinomio de mayor grado que la
fórmula integra exactamente. La Figura 4.21 es también de utilidad inmediata
4.40
SÓLIDOS BIDIMENSIONALES
27
α1 = 0.8168475730 ; β1 = 0.0915762135 ; γ1 = − ; 2γ3 = 0.1099517437
96
25
α2 = 0.1081030182 ; β2 = 0.4459484909 ; γ2 = ; 2γ4 = 0, 2233815897
96
4.41
Cálculo de Estructuras por el Método de Elementos Finitos
" "
" (e) "
"J " dξdη
A(e)
4.42
SÓLIDOS BIDIMENSIONALES
4.43
Cálculo de Estructuras por el Método de Elementos Finitos
np nq np nq (4.117)
" "
" " t
= BTi D Bj "J(e) "t Wṗ Wq̇ = " "
" (e) " Gij Wp Wq
p,q "J "
p=1 q=1 p=1 q=1
p,q
1 1−β
(e)
Kij = BTi DBj |J(e) |tdαdβ =
0 0
np np
T (e) t
= [Bi DBj |J |t]p Wp = [ (e) Gij ]pWp (4.118)
p=1 p=1 |J |
4.44
SÓLIDOS BIDIMENSIONALES
Figura 4.24 Viga en voladizo analizada con cuatro elementos Serendı́pitos de 8 nodos.
Valores del esfuerzo cortante en las secciones correspondientes a la
cuadratura 2 × 2 de Gauss-Legendre y extrapolación lineal a los nodos.
4.45
Cálculo de Estructuras por el Método de Elementos Finitos
4.46
SÓLIDOS BIDIMENSIONALES
(e)
Sustituyendo (4.123) en la expresión del vector fti de (4.120), se tiene
(e) τ cos β − σ sen β τ dx − σdy
fti = (e) Ni t ds = Ni t (4.124)
S σ cos β + τ sen β s(e) σdx + τ dy
Por otra parte, en el contorno en cuestión
∂x ∂y
dx = dξ = J11 dξ ; dy = dξ = J12 dξ (4.125)
∂ξ ∂ξ
+1 np
(e) τ J11 − σJ12 τ J11 − σJ12
fti = Ni t dξ = Ni t Wp (4.126)
−1 σJ11 + τ J12 σJ11 + τ J12
p=1 p
4.47
Cálculo de Estructuras por el Método de Elementos Finitos
con ella misma. Otros ejemplos de simetrı́a rotacional son las estructuras de
revolución y las estructuras con simetrı́a cı́clica que se estudian más adelante.
Una estructura simétrica puede tener cargas simétricas o antisimétricas. Un
sistema de cargas es antisimétrico si una reflexión de la estructura, con sus cargas,
seguida de un cambio de signo de todas ellas, produce la coincidencia con el
estado inicial (Figura 4.27c). En ambos casos, simétrico y antisimétrico, basta
con analizar la mitad simétrica de la estructura con las condiciones de contorno
siguientes:
Carga simétrica
Carga antisimétrica
4.48
SÓLIDOS BIDIMENSIONALES
♣ Ejemplo 4.2 Imponer las condiciones de contorno por simetrı́a en las mallas de
la Figura 4.28.
– Solución
Malla 1 : Viga en tensión plana (Figura 4.28a). Dadas las condiciones de
antisimetrı́a de la malla, basta con analizar la mitad de la misma con la
condición de desplazamiento vertical nulo en los nodos sobre el eje de
simetrı́a A − A .
Malla 2 : Placa bajo cargas puntuales (Figura 4.28b). La doble simetrı́a de malla
y cargas permite analizar un cuadrante de malla con las condiciones de
giros nulos sobre los lados 4-5 y 2-5 que se muestran.
4.49
Cálculo de Estructuras por el Método de Elementos Finitos
Figura 4.29 Elemento triangular de tres nodos apoyado sobre un terreno elástico
con coeficiente de balasto k.
4.50
SÓLIDOS BIDIMENSIONALES
donde todas las matrices y vectores tienen las expresiones usuales a excepción de
la matriz G(e) que se obtiene por
(e) 0 0
Gij = k(e) dx (4.131)
l(e) 0 Ni Nj
i j k
.. ..
0 0 . 0 0 . 0 0
i
.. ..
(e) 0 1 . 0 1/2 . 0 0
klij .. ..
G(e) = 0 0 . 0 0 . 0 0 (4.132)
3 .. .. j
0 1/2 . 0 1 . 0 0
.. ..
0 0 . 0 0 . 0 0
.. ..
0 0 . 0 0 . 0 0 k
4.51
Cálculo de Estructuras por el Método de Elementos Finitos
4.52
CAPÍTULO 5
SÓLIDOS DE REVOLUCIÓN
5.1 INTRODUCCIÓN
5.1
Cálculo de Estructuras por el Método de Elementos Finitos
5.2
SÓLIDOS DE REVOLUCIÓN
∂u ∂w ∂u ∂w
εr = ; εz = ; γrz = + (5.2)
∂r ∂z ∂z ∂r
2π(r + u) − 2πr u
εθ = = (5.3)
2πr r
5.3
Cálculo de Estructuras por el Método de Elementos Finitos
5.4
SÓLIDOS DE REVOLUCIÓN
σ = D (εε − ε0 ) + σ 0 (5.6)
5.5
Cálculo de Estructuras por el Método de Elementos Finitos
5.6
SÓLIDOS DE REVOLUCIÓN
donde
N1 0 | N2 0 | N3 0
N = [N1 , N2, N3 ] = (5.14)
0 N1 | 0 N2 | 0 N3
y
(e)
a1
T
a(e) = a
(e) = u1 , w1 , u2 , w2 , u3, w3 (5.15)
2
(e)
a3
∂Ni
0
∂r (e)
∂Ni 1
a
3 0
ε= ∂z
ui = B1 , B2 , B3 (e) = B a(e) (5.16)
w a 2
i=1
Ni
0
i
r (e)
a3
∂Ni ∂Ni
∂z ∂r
5.7
Cálculo de Estructuras por el Método de Elementos Finitos
3
σ= Bi D a(e) − D ε0 + σ 0 (5.19)
i=1
5.8
SÓLIDOS DE REVOLUCIÓN
(e)
Kij = 2π BTi D Bj r drdz (5.21)
A(e)
2×2 2×4 4×4 4×2
(d11 bi bj + d44 ci cj )r+ (d12 bi cj + d44 ci bj )r+
+2A(e) (d13 bi Nj + d31 bj Ni )+ +2A(e) d32 Ni cj
N N
(e) π +4(A(e) )2 d33 ir j
Kij =
drdz
2(A(e) )2 A(e)
(d21 ci bj + d44 bi cj )r+ (d22 cicj + d44 bi bj )r
+2A(e) d23 ci Nj
(e)
donde (·) indica valores en el baricentro del elemento. La expresión de Kij en este
caso se puede deducir directamente de la Figura 5.7 sustituyendo r por r, Ni por
1 , y el valor de la integral por el integrando multiplicado por el área del elemento.
3
(e)
Es fácil observar que esta integración es exacta para todos los términos de Kij a
NN
excepción del término ir j que se evalúa por defecto. No obstante, dicho error
no incide en el buen comportamiento del elemento, con el que pueden obtenerse
excelentes resultados utilizando una discretización suficientemente tupida en las
zonas donde se prevean mayores gradientes de tensiones [Z3],[Z8].
5.9
Cálculo de Estructuras por el Método de Elementos Finitos
f (e) = 2π (e)
NT b rdA + 2π (e)
NT t rds + 2π BT D ε0 rdA−
6×1 A l A(e)
(e) (e) (e) (e)
− 2π BT σ 0 r dA = fb + ft + fε + fσ (5.23)
A(e)
donde la primera integral corresponde al vector de fuerzas másicas de volumen
(e) (e)
fb ; la segunda al de fuerzas de superficie ft ; la tercera al de fuerzas debidas a
(e)
deformaciones iniciales fε ; y la cuarta al de fuerzas debidas a tensiones iniciales
(e)
fσ . De nuevo la ec. (5.23) es válida para cualquier tipo de elemento.
5.10
SÓLIDOS DE REVOLUCIÓN
(e)
(d11 + d12 + d13 )bi r + 2A3 (d31 + d32 + d33 )
(e)
fεi = π α∆T
(5.27)
(d21 + d22 + d23 )ci r
siendo r la coordenada nodal del baricentro del elemento y dij los elementos de la
matriz constitutiva.
(e)
La expresión de fσi para el elemento triangular de tres nodos sometido a
tensiones iniciales constantes se obtiene de forma exacta por
(bi σr0 + ci τrz
0
)r + 23 A(e) σθ0
(e)
fσi = −π (5.29)
0 0
(ci σz + bi τrz )r
con
(e) Ri
pi = (5.30b)
Zi
5.11
Cálculo de Estructuras por el Método de Elementos Finitos
con
Ni 0 (e) ri
Ni = ; xi = (5.32)
0 Ni zi
+1 +1 n
(e)
Kij = 2π BTi DBj r dr dz = 2π BTi DBj Nk rk J(e) dξ dη
A(e) −1 −1
k=1
(5.33)
y el vector de fuerzas nodales equivalentes másicas
+1 +1 n
Nk rk J(e)
(e)
fb = 2π NTi br dr dz = 2π T
Ni b dξ dη (5.34)
i A(e) −1 −1 k=1
5.12
SÓLIDOS DE REVOLUCIÓN
5.5 CONCLUSIONES
5.13
CAPÍTULO 6
SÓLIDOS TRIDIMENSIONALES
6.1 INTRODUCCIÓN
6.1
Cálculo de Estructuras por el Método de Elementos Finitos
donde u, v, w son los desplazamientos del punto según los ejes cartesianos x, y, z,
respectivamente.
6.2
SÓLIDOS TRIDIMENSIONALES
con
∂u ∂v ∂w
εx = ; εy = ; εz =
∂x ∂y ∂z
(6.3)
∂u ∂v ∂u ∂w ∂v ∂w
γxy = + ; γxz = + ; γyz = +
∂y ∂x ∂z ∂x ∂z ∂y
donde εx , εy , εz son las deformaciones normales y γxy , γxz , γyz las deformaciones
tangenciales.
T
σ = σx, σy , σz , τxy , τxz , τyz (6.4)
donde σx , σy , σz son las tensiones normales y τxy , τxz ,τyz son las tensiones
tangenciales. En la Figura 6.3 se muestra el convenio de signos de dichas tensiones.
6.3
Cálculo de Estructuras por el Método de Elementos Finitos
σ = D (εε − ε0 ) + σ 0 (6.5)
donde V y A son el volumen y la superficie del cuerpo sobre los que actúan las
fuerzas de masa b, de superficie t y puntuales qi , respectivamente.
Como en los problemas de elasticidad estudiados en los dos capı́tulos
anteriores, en la expresión del PTV intervienen sólo primeras derivadas de
los desplazamientos, lo que exige simplemente continuidad de clase C0 a la
aproximación de elementos finitos.
6.4
SÓLIDOS TRIDIMENSIONALES
u
N1 u1 + N2 u2 + N3 u3 + N4 u4
4
(e)
u= v = N1 v1 + N2 v2 + N3 v3 + N4 v4 = Ni ai = N a(e)
w N1 w1 + N2 w2 + N3 w3 + N4 w4 i=1
(6.9)
donde
y
(e)
a 1
(e)
a2
a(e) = (6.12)
(e)
a
3
(e)
a4
ui
(e)
ai = vi (6.13)
wi
6.5
Cálculo de Estructuras por el Método de Elementos Finitos
u = α1 + α2 x + α3 y + α4 z
v = α5 + α6 x + α7 y + α8 z (6.14)
w = α9 + α10 x + α11 y + α12 z
Las constantes αi se obtienen sustituyendo adecuadamente las coordenadas
de los nodos e igualando los desplazamientos a sus valores nodales. Como
hemos utilizado la misma aproximación para todos los desplazamientos, basta
con calcular las cuatro constantes para un solo desplazamiento. Ası́, considerando
el desplazamiento u
u1 = α1 + α2 x1 + α3 y1 + α4 z1
u2 = α1 + α2 x2 + α3 y2 + α4 z2
(6.15)
u3 = α1 + α2 x3 + α3 y3 + α4 z3
u4 = α1 + α2 x4 + α3 y4 + α4 z4
6.6
SÓLIDOS TRIDIMENSIONALES
de donde se deduce, por comparación con (6.9), que la función de forma del nodo
i es
1
Ni = (ai + bi x + ci y + di z) (6.17)
6V (e)
6.7
Cálculo de Estructuras por el Método de Elementos Finitos
B = [B1, B2 , B3 , . . . , Bn ] (6.20)
B = [B1 , B2 , B3 , B4 ] (6.22)
donde K(e) es la matriz de rigidez del elemento, f (e) es el vector de fuerzas nodales
equivalentes y q(e) es el vector de fuerzas nodales de equilibrio (que desaparece en
el ensamblaje). La matriz de rigidez tiene la expresión usual
K(e) = BT D B dV (6.25)
V (e)
3n × 3n 3n × 6 6 × 6 6 × 3n
6.8
SÓLIDOS TRIDIMENSIONALES
(e)
La expresión desarrollada de Kij para este elemento se presenta en la
Figura 6.5.
6.9
Cálculo de Estructuras por el Método de Elementos Finitos
Es fácil encontrar una forma explı́cita del vector de fuerzas nodales equivalentes
del elemento tetraédrico de cuatro nodos como a continuación se muestra.
Fuerzas de volumen
(e)
b1
f
(e)
fb
(e) 2
fb = = NT b dV (6.29)
(e) V (e)
fb
3
f (e)
b4
con
(e) bx
(e) V
fb = NTi b dV = by (6.30)
i V (e) 4
bz
Fuerzas de superficie
ft1
(e)
t2
f
(e)
ft = = NT tdA (6.31)
ft3 A(e)
(e)
(e)
ft4
con
(e)
fti = NTi t dA (6.32)
A(e)
6.10
SÓLIDOS TRIDIMENSIONALES
(e)
(e) A
ft = 234 [0, 0, 0, tx, ty , tz , tx , ty , tz , tx , ty , tz ]T (6.35)
3
y
Fuerza actuando sobre la cara definida por los nodos 1-3-4
(e)
(e) A
ft = 134 [tx , ty , tz , 0, 0, 0, tx, ty , tz , tx , ty , tz ]T (6.36)
3
(e)
f ε1
(e)
fε2
(e)
fε = = BT D ε0 dV (6.37)
(e)
V (e)
f ε3
(e)
fε4
donde
0 0
(d11 εx + d12 εy + d13 ε0z )bi
(e) 1
fεi = BTi D ε0 dV = (d21 ε0x + d22 ε0y + d23 ε0z )ci (6.38)
V (e) 6
(d31 ε0x + d32 ε0y + d33 ε0z )di
(e) E (1 − ν)α(∆T )
fεi = [b , c , d ]T (6.39)
6(1 + ν)(1 − 2ν) i i i
6.11
Cálculo de Estructuras por el Método de Elementos Finitos
donde
bi σx0 + ci τxy
0 0
+ di τxz
(e)
fσi = − B σ 0 dV = − 1 ci σy0 + bi τxy
0 0
+ di τyz (6.41)
i
V (e) 6
di σz0 + bi τxz
0 0
+ ci τyz
(x − xc ) (y − yc ) (z − zc )
ξ= ; η= ; ζ= (6.42)
a b c
donde (xc , yc , zc ) son las coordenadas del centro de gravedad del elemento. De
(6.42) se deduce
dξ 1 dη 1 dζ 1
= ; = ; = (6.43)
dx a dy b dz c
6.12
SÓLIDOS TRIDIMENSIONALES
Por ser el elemento recto, las derivadas cartesianas de las funciones de forma
se pueden calcular directamente por la expresión
Ni (ξ, η, ζ) = 1 (6.47b)
i=1
6.13
Cálculo de Estructuras por el Método de Elementos Finitos
donde lIi (ξ) es el polinomio de Lagrange de grado I en ξ que pasa por el nodo
i, etc. Como en el caso de elementos rectangulares, es usual escoger la misma
aproximación polinómica de Lagrange en cada una de las tres direcciones ξ, η y ζ.
Los términos polinómicos de las funciones de forma se obtienen sencillamente
del tetraedro de Pascal en tres dimensiones. En la Figura 6.7 se muestran los
elementos de 8 y 27 nodos de esta familia y los términos polinómicos del tetraedro
de Pascal que intervienen en sus funciones de forma cuya deducción se detalla en
los dos subapartados siguientes.
1
Ni (ξ, η, ζ) = (1 + ξi ξ) (1 + ηi η) (1 + ζi ζ) (6.49)
8
Adviértase que:
6.14
SÓLIDOS TRIDIMENSIONALES
por elemento. Debido a ello, suele ser usual su utilización incorporando algunas
modificaciones para mejorar su comportamiento a flexión. En la referencia [O3] se
describen algunas de estas técnicas.
6.15
Cálculo de Estructuras por el Método de Elementos Finitos
Nodos laterales
1 1
Ni = ηi2 (η 2 − ηηi )ζ 2 (ζ 2 − ζζi )(1 − ξ 2 ) + ζi2 (ζ 2 − ζζi )+ 2, 4, 6, 8
4 4 ; i = 10, 12, 14, 16
2 2 2 1 2 2 2 2 2
+ ξi (ξ − ξξi )(1 − η ) + ξi (ξ − ξξi)ηi (η − ηηi )(1 − ζ ) 20, 22, 24, 26
4
(6.51)
6.16
SÓLIDOS TRIDIMENSIONALES
6.17
Cálculo de Estructuras por el Método de Elementos Finitos
6.18
SÓLIDOS TRIDIMENSIONALES
1 1, 3, 5, 7
Ni = (1 + ξi ξ)(1 + ηi η)(1 + ζi ζ)(ξi ξ + ηi η + ζi ζ − 2) ; i = (6.54)
8 13, 15, 17, 19
6.19
Cálculo de Estructuras por el Método de Elementos Finitos
Nodos laterales
1
Ni = (1 − ξ 2 )(1 + ηi η)(1 + ζi ζ) ; i = 2, 6, 14, 18
4
1
= (1 − η2 )(1 + ζi ζ)(1 + ξi ξ) ; i = 4, 8, 16, 20 (6.55)
4
1
= (1 − ζ 2 )(1 + ηi η)(1 + ξi ξ) ; i = 9, 10, 11, 12
4
Adviértase que:
1 Todas las funciones de forma contienen un polinomio completo de segundo
grado más los términos ξη2 , ξ 2 η, ξ 2 ζ, ξ.ζ 2 , ζ 2 η, η2 ζ, ξηζ, ξ 2 ηζ, ξη2 ζ y ξηζ 2 .
2 Las funciones de forma satisfacen las condiciones (6.47).
Conviene destacar que el elemento hexaédrico Serendı́pito de 20 nodos tiene la
misma aproximación cuadrática que el Lagrangiano de 27 nodos con siete nodos
menos, lo que representa un ahorro total de 21 variables nodales. Esto explica la
mayor popularidad del primero para análisis tridimensionales. Estas diferencias
son incluso más acusadas para elementos de órdenes superiores [O3].
L1 + L2 + L3 + L4 = 1 (6.57)
x= Li x i , y= Li yi , z = Li zi (6.58)
i=1 i=1 i=1
6.20
SÓLIDOS TRIDIMENSIONALES
Estas tres ecuaciones junto con la (6.57) permiten eliminar las Li en función
de las coordenadas cartesianas. Es fácil comprobar que
l
Li = (ai + bi x + ci y + di z) = Ni (6.59)
6V (e)
6.21
Cálculo de Estructuras por el Método de Elementos Finitos
donde los coeficientes ai, bi , ci , di coinciden con los de la ec.(6.17) para las funciones
de forma del elemento tetraédrico de cuatro nodos. Se deduce, por tanto, que las
coordenadas de volumen coinciden con las funciones de forma de dicho elemento.
Obsérvese que de (6.59) se pueden obtener las derivadas de las coordenadas de
volumen con respecto a las cartesianas como
x − x1 y − y1 z − z1
α= ; β= ; γ= (6.61)
a b c
dα 1 dβ 1 dγ 1
= ; = ; = (6.62)
dx a dy b dz c
dx dy dz = abc dα dβ dγ (6.63)
1 1−α 1−β−γ
(e)
f (x, y, z)dx dy dz = f (α, β, γ)abc dα dβ dγ (6.64)
V 0 0 0
6.22
SÓLIDOS TRIDIMENSIONALES
Es fácil encontrar que las funciones de forma del elemento tetraédrico de cuatro
nodos se expresan en coordenadas naturales como
N1 = 1 − α − β − γ ; N2 = α; N3 = β; N4 = γ (6.65)
donde los subı́ndices I, J, K, y L coinciden con los exponentes que afectan a cada
coordenada de volumen en la expresión de la función de forma Ni , cumpliéndose,
obviamente, que I + J + K + L = M, siendo M el grado del mayor polinomio
completo contenido en Ni . Por otra parte, lIi (Lj ) es el polinomio de Lagrange
de grado I en la variable Lj que pasa por el nodo i (ver ec.(3.6)). Dadas las
caracterı́sticas tridimensionales del tetraedro la asignación de las coordenadas
I, J, K, L a cada nodo es algo más complicada que en el caso bidimensional, como
puede apreciarse en la Figura 6.15.
6.23
Cálculo de Estructuras por el Método de Elementos Finitos
Nodo 1
Posición (I, J, K, L) : (2, 0, 0, 0). Coordenadas de volumen : (1, 0, 0, 0)
L1 − 12 L1
N1 = l21 (L1 ) = = (2L1 − 1)L1 (6.68)
1 − 12
Nodo 2
L L
N2 = l12 (L1 ) l12 (L2 ) = 11 12 = 4L1 L2 (6.69)
2 2
N3 = (2L2 − 1) ; N7 = 4 L2 L4
N4 = 4L2 L3 ; N8 = 4 L3 L4
(6.70)
N5 = (2 L3 − 1)L − 3 ; N9 = 4 L1 L4
N6 = 4 L1 L3 ; N10 = (2 L4 − 1) L4
6.24
SÓLIDOS TRIDIMENSIONALES
x
n
xi
x= y = Ni yi = N x(e) (6.71)
z i=1 zi
6.25
Cálculo de Estructuras por el Método de Elementos Finitos
con
Ni
N = [N1 , N2, . . . , Nn ] ; Ni =
Ni
; Ni = f (ξ, η, ζ) (6.72)
Ni
donde Ni es la misma función de forma utilizada para interpolar el campo de
desplazamientos.
La ec.(6.71) expresa una relación entre las coordenadas cartesianas y las
naturales. Dicha relación es biunı́voca si se cumple que el determinante del
6.26
SÓLIDOS TRIDIMENSIONALES
donde J(e) es la matriz Jacobiano cuyos términos se calculan haciendo uso de las
relaciones isoparamétricas (6.72) en la forma
∂Ni ∂Ni ∂Ni
∂ξ xi ∂ξ yi z
∂ξ i
n
∂Ni ∂Ni ∂Ni
J(e) = ∂η xi ∂η yi ∂η i
z (6.74)
i=1
∂Ni ∂Ni ∂Ni
∂ζ xi ∂ζ yi z
∂ζ i
6.27
Cálculo de Estructuras por el Método de Elementos Finitos
donde
(e)
J 1k
b̄i
3
(e) ∂Ni
c̄i = (6.78)
J 2k
k=1
∂ξk
d¯i
(e)
J 3k
ξ1 = ξ, ξ2 = η y ξ3 = ζ.
La expresión genérica de la matriz de rigidez de un elemento tridimensional en
coordenadas naturales es, por tanto,
+1 +1 +1
(e)
Kij = BTi D Bj dV = BTi (ξ, η, ζ)D Bj (ξ, η, ζ)J(e) dξ dη dζ =
V (e) −1 −1 −1
+1 +1 +1
= Gij (ξ, η, ζ) dξ dη dζ (6.79a)
−1 −1 −1
con
(d11 b̄ij + d44 c̄ij + d55 d¯ij ) (d12 b̄i c̄j + d44 c̄i b̄j ) (d13 b̄i d¯j + d55 d¯i b̄j )
Gij = (d21 c̄i b̄j + d44 b̄i c̄j ) (d22 c̄ij + d44 b̄ij + d66 d¯ij ) (d23 c̄i d¯j + d66 d¯i c̄j ) J(e)
(d31 d¯i b̄j + d55 b̄i d¯j ) (d32 d¯i c̄j + d66 c̄i d¯j ) (d33 d¯ij + d55 b̄ij + d66 c̄ij )
(6.79b)
y b̄ij = b̄i b̄j , c̄ij = c̄i c̄j y d¯ij = d¯i d¯j , donde b̄i , c̄i , d¯i se han definido en (6.78) y dij
son las componentes de la matriz D de (6.6). De las expresiones de la inversa del
Jacobiano se deduce que en la matriz G intervienen expresiones racionales [O3],
por lo que su integración analı́tica es sumamente complicada y hay que recurrir a
la integración numérica.
En elementos tetraédricos que utilicen coordenadas de volumen la interpolación
isoparamétrica se define similarmente a la ec.(6.71), siendo ahora Ni =
f (L1 , L2 , L3 , L4). Si los tetraedros son de lados rectos el cálculo de las derivadas
cartesianas de las funciones de forma es inmediato y las integrales sobre el elemento
pueden calcularse exactamente [O3]. Si el elemento es de lados curvos es más
conveniente, como ocurrı́a en el caso bidimensional, operar con las coordenadas
naturales, lo que simplemente implica sustituir L2, L3 y L4 por α, β y γ,
respectivamente, y L1 por 1-α-β-γ. A partir de aquı́ el cálculo de las componentes
de la matriz Bi sigue idénticos pasos a los explicados entre las ecs.(6.75)-(6.79)
para elementos hexaédricos, sin más que sustituir las coordenadas ξ, η, ζ por α, β, γ,
respectivamente.
Por consiguiente, la matriz de rigidez del elemento tetraédrico isoparamétrico
de lados curvos tiene una expresión similar a la (6.79a)
1 1−α 1−α−β
(e)
Kij = Gij (α, β, γ)dα dβ dγ (6.80)
0 0 0
6.28
SÓLIDOS TRIDIMENSIONALES
6.29
Cálculo de Estructuras por el Método de Elementos Finitos
Todas las consideraciones sobre selección del orden de integración hechas para
elementos bidimensionales siguen siendo válidas para el caso tridimensional y no
merecen mayor comentario.
+1 +1 +1
(e)
Kij = BTi D Bj dx dy dz = BTi D Bj J(e) dξ dη dζ =
V (e) −1 −1 −1
np nq nr
(e)
np nq n r
T
= Bi D Bj J Wp Wq Wr = [Gij ]p,q,r Wp Wq Wr (6.85)
p=1 q=1 r=1 p,q,r p=1 q=1 r=1
+1 +1 +1
(e)
fi = (e)
NTi b dx dy dz = NTi b J(e) dξ dη dζ =
V −1 −1 −1
np nq nr
= NTi bJ(e) Wp Wq Wr (6.86)
p=1 q=1 r=1 p,q,r
6.30
SÓLIDOS TRIDIMENSIONALES
2 ; δ= 3
α = 0.58541020 ; β = 0.13819660 ; γ = − 15 40
El caso de fuerzas sobre una de las caras del elemento es algo más complicado.
Para explicar el proceso consideremos que actúa una fuerza tn ortogonal a la cara
situada en ζ = +1 y definida por los nodos 5 a 8 (ver Figura 6.19). Para el
cálculo del vector de fuerzas de superficie precisamos conocer el término t dA,
donde dA es el diferencial de área en dicha cara, y t el vector de fuerzas en ejes
globales actuantes sobre la superficie en cuestión. Ası́, si nx , ny , nz son los cosenos
directores de la normal a la superficie, se cumple
6.31
Cálculo de Estructuras por el Método de Elementos Finitos
31 × V
V 32
n= (6.89)
|V
31 × V
32 |
31 × V
y recordando que dA = |V 32 |, se tiene que
(e)
J12 J23 − J22 J13
1 1 (e)
n= J21 J13 − J11 J23 dξ dη = j dξ dη (6.90)
dA
dA
J11 J32 − J21 J12 ζ=+1
(e)
donde las Jij se obtienen de (6.74).
Por consiguiente, la expresión final del vector de fuerzas de superficie es
+1 +1
(e)
fti = Ni tn n dA = Ni tn j(e) dξ dη =
A(e) −1 −1
np nq
6.32
SÓLIDOS TRIDIMENSIONALES
Figura 6.20 Análisis de la flexión de una viga con diferentes elementos de sólido
tridimensionales.
6.33
Cálculo de Estructuras por el Método de Elementos Finitos
Por otra parte, hay que señalar que la discretización de un sólido tridimensional
en tetraedros es mucho más compleja que en el caso de los más sencillos triángulos
en dos dimensiones. El problema se simplifica si se parte de una discretización en
hexaedros y se divide cada uno de ellos en cinco elementos tetraédricos [P1]. No
obstante, dada la superioridad de los elementos hexaédricos antes mencionada,
este procedimiento no parece apropiado desde el punto de vista de obtener
una mayor precisión en el cálculo. Hoy en dı́a se ha avanzado mucho en las
técnicas de generación de mallas tetraédricas, lo que ha favorecido la utilización de
estos elementos en el análisis de geometrı́as tridimensionales complejas utilizando
técnicas de remallado adaptable [B7], [G5], [O3], [P2], [P3], [P4].
6.34
CAPÍTULO 7
FLEXIÓN DE VIGAS
7.1 INTRODUCCIÓN
7.1
Cálculo de Estructuras por el Método de Elementos Finitos
u(x, y, z) = − zθ(x)
v(x, y, z) = 0 (7.1)
w(x, y, z) = w(x)
7.2
FLEXIÓN DE VIGAS
dw dw
θ = y u = −z (7.2)
dx dx
du d2 w
εx = = −z
dx dx2 (7.3)
εy = εz = γxy = γxz = γyz = 0
d2 w
σx = E εx = − z E (7.4)
dx2
d2w d2 w
M = − zσx dA = z2 E dA = EI = EIχ (7.5)
A A dx2 dx2
7.3
Cálculo de Estructuras por el Método de Elementos Finitos
w = αo + α1 x + α2 x2 + α3 x3 (7.8)
w1 = αo + α1 x1 + α2 x21 + α3 x31
dw
= α1 + 2α2x1 + 3α3 x21
dx 1 (7.9)
w2 = αo + α1 x2 + α2 x22 + α3 x32
dw
= α1 + 2α2x2 + 3α3 x22
dx 2
Una vez resuelto el sistema anterior se puede reescribir (7.8), tras sustituir
convenientemente las expresiones de las αi , como
l(e) dw l(e) dw
w = N1 w1 + N 1 ( )1 + N2 w2 + N 2 ( ) (7.10)
2 dx 2 dx 2
7.4
FLEXIÓN DE VIGAS
1 1
N1 = (2 − 3ξ + ξ 3 ) ; N2 = (2 + 3ξ − ξ 3 )
4 4 (7.11)
1 1
N 1 = (1 − ξ − ξ 2 + ξ 3 ) ; N 2 = (−1 − ξ + ξ 2 + ξ 3 )
4 4
w = N a(e) (7.13)
donde
dw dw T
N = N1 , N 1 , N2 , N 2 y a(e) = w1 , , w2 , (7.14)
dx 1 dx 2
son la matriz de funciones de forma y el vector de movimientos (desplazamientos
y giros) nodales del elemento, respectivamente.
7.5
Cálculo de Estructuras por el Método de Elementos Finitos
l(e) dw 2 dw d2 w 4 d2 w
dx = dξ ; = (e) y = (7.15)
2 dx l dξ dx2 2
(l(e) ) dξ
2
d2 w 4 d2 N1 l(e) d2 N 1 dw d2 N2 l(e) d2 N 2 dw
χ = 2 = (e) w1 + + w2 + =
dx (l )2 dξ 2 2 dξ 2 dx 1 dξ 2 2 dξ 2 dx 2
w1
6ξ (−1 + 3ξ) −6ξ (1 + 3ξ) dw
= , , ,
dx 1 = Bf a(e) (7.16)
(l(e) )2 l(e) (l(e) )2 l(e) w
dw2
dx 2
+1 T l(e)
δχ EIχdx = δa(e) BTf (EI) Bf dξ a(e) =
l( e) −1 2
+1 T 2 2 dw
(7.17)
ql(e)
=− δa(e) NT dξ + δwi Zi + δ M
−1 2 i=1 j=1 dx j j
7.6
FLEXIÓN DE VIGAS
donde la matriz de rigidez del elemento de viga puede calcularse de forma explicita
por
12 6l(e) −12 6l(e)
+1 (e)
.. 2 2
EIl (e) EI . 4(l(e) ) −6l(e) 2(l(e) )
K(e) = BT B dξ =
...
−1 2 l3 12 −6l(e)
... (e) 2
sim. 4(l )
(7.19)
El lector reconocerá en las componentes del vector f (e) de (7.20) los valores,
con los signos de la Figura 7.4, de las reacciones verticales y los momentos en los
extremos de una viga biempotrada bajo carga uniforme. Esta coincidencia es, no
obstante, un caso muy particular, debido a las caracterı́sticas especiales de la carga
uniforme, no siendo por tanto extrapolable a otro tipo de cargas ni de elementos
[O3].
Una vez obtenidos los desplazamientos y los giros nodales se puede obtener el
momento flector en cualquier punto del elemento por la expresión
M = EI χ = EI B a(e) (7.22)
7.7
Cálculo de Estructuras por el Método de Elementos Finitos
7.8
FLEXIÓN DE VIGAS
dw du dw
γxz = + = −θ = − φ (7.25)
dx dz dx
Por consiguiente, 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 dos tensiones no nulas σx y τxz se relacionan con las correspondientes
deformaciones por
dθ
σx = Eεx = − zE = −zEχ
dx (7.26)
dw
τxz = Gγxz = G −θ
dx
dθ .
donde G es el módulo de rigidez y χ = dx
El momento flector y el esfuerzo cortante se definen, de acuerdo con los signos
de la Figura 7.6, como
dθ
M =− zσx dA = EI = EIχ
A dx (7.27)
dw
Q= τxz dA = GA − θ = GAγxz
A dx
Obsérvese que la variación de σx con el canto es lineal, lo cual puede
considerarse como “exacto” dentro de la hipótesis de la teorı́a de vigas. Por el
7.9
Cálculo de Estructuras por el Método de Elementos Finitos
y
Q = α A G γxz = A∗ G γxz (7.29)
l p
q
(δεx σx + δγxz τxz )dV =− δwqdx − δwi Zi + δθj Mj (7.30)
V 0 i=1 j=1
7.10
FLEXIÓN DE VIGAS
Es fácil ver que haciendo uso de las ecs.(7.24) - (7.29) el primer miembro de
(7.30) puede modificarse como
dθ dw
−zσx δ + τxz δ −θ dV =
V dx dx
l
= δχ −zσx dA + δγxz τxz dA dx =
0 A A
l (7.31)
= δχM + δγxz Q dx =
0
l
dθ
dθ dw dw
= δ EI +δ − θ GA∗ − θ dx
0 dx dx dx dx
7.11
Cálculo de Estructuras por el Método de Elementos Finitos
χ = Bf a(e)
(7.35)
γxz = Bc a(e)
donde
2 dN1 2 dN2 1 1
Bf = 0, (e) , 0, (e) = 0, − (e) , 0 (e)
l dξ l dξ l l
2 dN1 2 dN2 1 −(1 − ξ) 1 −(1 + ξ)
Bc = , −N1 , , −N2 = − , , ,
l(e) dξ l(e) dξ l(e) 2 l(e) 2
(7.36)
son las matrices de deformación de flexión y cortante del elemento, y
a(e) = [w1 , θ1, w2 , θ2 ]T (7.37)
es el vector de movimientos nodales del elemento.
La expresión de los trabajos virtuales (7.30) puede escribirse, haciendo uso de
las ecs.(7.31) - (7.37), como
T
δa (e) T T ∗
Bf (EI)Bf + Bc (GA )Bc dx a(e) =
l ( e)
T
(7.38)
= δa(e) N
T
(−q)dx + δa (e) T q(e)
l(e)
7.12
FLEXIÓN DE VIGAS
donde
(e) (e)
K(e) = Kf + Kc (7.40)
y
BTc (GA∗ )Bc dx
(e) (e)
Kf = BTf (EI)Bf dx ; Kc = (7.41)
l(e) l(e)
son las matrices de rigidez correspondientes a los efectos de flexión y cortante cuya
suma es la matriz de rigidez total del elemento;
T
f (e) =− N q dx, con N = [N1 , 0, N2 , 0] (7.42)
l(e)
el vector de fuerzas nodales equivalentes debidas a la carga repartida q; y
q(e) = [P1, M1 , P2 , M2 ]T (7.43)
el vector de fuerzas nodales de equilibrio que permite ensamblar las contribuciones
de los distintos elementos en la matriz de rigidez y en el vector de fuerzas globales.
Todas las integrales anteriores pueden transformarse sobre el dominio
normalizado del elemento.
(e)
Ası́, teniendo en cuenta que dx = l 2 dξ, las ecs.(7.41) y (7.42) se escriben como
+1 +1
l(e) l(e)
BTc (GA∗ ) Bc
(e) (e)
Kf = BTf (EI) Bf dξ ; Kc = dξ (7.44)
−1 2 −1 2
y
+1
T l(e)
f (e) = − N q dξ (7.45)
−1 2
Las integrales anteriores pueden evaluarse numéricamente por una cuadratura
unidimensional de Gauss-Legendre (Apartado 3.4).
Adviértase que la matriz de rigidez del elemento puede también obtenerse por
la expresión general
K(e) = BT D B dx (7.46)
l(e)
donde !
Bf EI 0
B= y D= (7.47)
Bc 0 GA∗
7.13
Cálculo de Estructuras por el Método de Elementos Finitos
GA∗ GA∗ ∗
GA∗
l 2 − GA
l 2
GA∗ EI ∗
GA∗
3 l + − GA 6 l − EI
w1 V1
w1 = 0
l 2 l
θ1 M1 θ1 = 0
= (7.51)
.. GA∗
− GA
∗ w2 P
. l 2
.. θ2 0
. GA∗ EI
Simetr. 3
l + l
7.14
FLEXIÓN DE VIGAS
γ l l3
w2 = + P (7.54)
γ + 1 GA∗ 3EI
3
En el caso de una sección rectangular I = bh ∗ 5
12 , A = 6 bh y con ν = 0.25
h 2 3
γ=3 = 2 (7.55)
l λ
Es conocido que en 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 (7.57).
7.15
Cálculo de Estructuras por el Método de Elementos Finitos
l + l3 P
w2 γ GA∗ 3EI 3(4λ2 + 3)
ϕ= = 3 = (7.58)
f
(w2 )exacta γ +1 l 4λ2 (λ2 + 3)
3EI P
1 l(e) −1 l(e)
2 2 2 2
... l(e) (e) l(e)
GA∗ (e) −l 2
(e) 4 4
Kc = (7.59)
l .. l (e)
. 1 − 2
2
.. l(e)
Simetr. .
4
GA∗ − GA
∗ l l3
+ 4EI l2
K= l ∗ 2 y F= GA∗ 2EI (7.60)
− GA GA∗ l + EI l2 l
2 4 l 2EI EI
7.16
FLEXIÓN DE VIGAS
7.17
Cálculo de Estructuras por el Método de Elementos Finitos
7.18
FLEXIÓN DE VIGAS
l2 l2 GA∗ l3
Kf + Kc a = f = f̄ (7.64)
3 3EI 3EI
β Kc a = f (7.66)
2
donde β = 4Gα l
E h . En el lı́mite, para vigas infinitamente esbeltas h → 0 y
β→∞y
1
Kc a = f → 0 (7.67)
β
7.19
Cálculo de Estructuras por el Método de Elementos Finitos
2−1×2 =0
2−1×1 =1> 0
2−1×2 =0
7.20
FLEXIÓN DE VIGAS
7.4 CONCLUSIONES
7.21
CAPÍTULO 8
8.1 INTRODUCCIÓN
8.1
Cálculo de Estructuras por el Método de Elementos Finitos
u=v=0 (8.1)
En otras palabras, los puntos del plano medio sólo se mueven verticalmente.
2) Todos los puntos contenidos en una normal al plano medio tienen
aproximadamente el mismo desplazamiento vertical.
3) La tensión normal σz es despreciable.
4) Los puntos sobre rectas normales al plano medio antes de la deformación,
permanecen sobre rectas también ortogonales a la deformada del plano medio
después de la deformación.
Las hipótesis 1, 2 y 4, permiten definir el campo de desplazamientos a través del
espesor de la placa. La tercera hipótesis afecta a la relación tensión-deformación,
como se verá en el Apartado 8.2.4.
Figura 8.1 Definición geométrica de una placa y convenio de signos para despla-
zamientos y giros.
8.2
PLACAS DELGADAS Y GRUESAS
u = [w, θx , θy ]T (8.3)
Figura 8.2 Deformación del plano medio de una placa delgada y giro de la normal.
8.3
Cálculo de Estructuras por el Método de Elementos Finitos
∂u ∂ 2w
εx = = −z 2
∂x ∂x
∂v ∂ 2w
εy = = −z 2 ; εz 0
∂y ∂y
∂u ∂v ∂ 2w
γxy = + = −2z
∂y ∂x ∂x∂y
(8.6)
∂w ∂u ∂w ∂v
γxz = + = 0 ; γyz = + =0
∂x ∂z ∂y ∂z
8.4
PLACAS DELGADAS Y GRUESAS
donde
t3
D̂f = D (8.12)
12
es la matriz constitutiva de flexión y
T
∂ 2w ∂ 2w ∂ 2w
ε̂εf = − 2 , − 2 , −2 (8.13)
∂x ∂y ∂x∂y
8.5
Cálculo de Estructuras por el Método de Elementos Finitos
εf = zε̂εf (8.14)
+ 2t
= δε̂εTf zσ
σ dz dA = δε̂εTf σ
σ̂ f dA (8.16)
A − 2t A
8.6
PLACAS DELGADAS Y GRUESAS
∂Qx ∂Qy
+ +q =0 (8.19)
∂x ∂y
Equilibrio de momentos
∂My ∂Mxy
My = 0 ⇒ dy dx + dx dy−
∂y ∂x
∂Qy ∂Qx dy dy
− Qy + dy dx − dx dy + q dxdy =0
∂y ∂x 2 2
∂Mx ∂Mxy
Mx = 0 ⇒ dx dy + dy dx−
∂x ∂y
∂Qx ∂Qy dx dy
− Qx + dx dy − dy dx + q dxdy =0 (8.20)
∂x ∂y 2 2
8.7
Cálculo de Estructuras por el Método de Elementos Finitos
∂Mx ∂Mxy
+ − Qx = 0 (8.22)
∂x ∂y
∂ 2Mx ∂ 2Mxy ∂ 2 My
+ 2 + = −q (8.23)
∂x2 ∂x∂y ∂y
∂ 4w ∂ 4w ∂ 4w q
4 + 2 2 2 + 4 =
∂x ∂x ∂y ∂y D
q Et 3
o 4 w = D con D = 12(1−ν 2) (8.24)
que es una ecuación diferencial de cuarto orden que relaciona la flecha con la carga
repartida y las propiedades del material. Dicha ecuación con sus correspondientes
condiciones de contorno, es el punto de partida para resolver analı́ticamente
problemas de placas isótropas.
Una vez calculada la flecha por integración de (8.24) los momentos flectores en
cada punto se obtienen por (8.11). Por otra parte, sustituyendo (8.11) en (8.21) y
(8.22) se deducen las expresiones de los esfuerzos cortantes
3 3
∂ w ∂ 3w ∂ w ∂ 3w
Qx = −D + ; Qy = −D + (8.25)
∂x3 ∂x∂y2 ∂y3 ∂y∂x2
3 Qx 3 Qy
(τxz )max = ; τyz = (8.26)
2 t max 2 t
8.8
PLACAS DELGADAS Y GRUESAS
w = α1 + α2 x + α3 y + α4 x2 + α5 xy+....
(8.27)
(hasta 3n términos)
8.9
Cálculo de Estructuras por el Método de Elementos Finitos
w = α1 + α2 x + α3 y + α4 x2 + α5 xy + α6 y2 + α7 x3 + α8 x2 y+
(8.29)
+α9 xy2 + α10 y3 + α11 x3 y + α12 xy3
8.10
PLACAS DELGADAS Y GRUESAS
8.11
Cálculo de Estructuras por el Método de Elementos Finitos
con
2 2 2 ¯
∂ Ni
− ∂x2 − ∂∂xN̄2i − ∂∂xN̄2 i
2 2 2¯
B = [B1 , B2 , B3 , B4 ] , Bi = − ∂∂yN2i − ∂∂yN̄2i − ∂∂yN̄2 i
2 2 2 N̄
¯
∂ Ni ∂ N̄i
−2 ∂x∂y −2 ∂x∂y −2 ∂∂x∂y i
∂2 1 ∂2 ∂2 1 ∂2 ∂2 1 ∂2
= 2 2 ; = 2 y = (8.33)
∂x2 a ∂ξ ∂y2 b ∂η2 ∂x∂y ab ∂ξ∂η
Dado que las funciones de forma tienen una variación cúbica se deduce que las
curvaturas, y por consiguiente los momentos flectores varı́an linealmente dentro
del elemento MZC.
Siguiendo el procedimiento usual se tiene
8.12
PLACAS DELGADAS Y GRUESAS
(e)
Pi
Ni
(e) T
fi = Mxi = N i qdxdy = q N̄i dxdy (8.37)
A(e) A(e) ¯
Myi N̄ i
Por otra parte, q(e) en (8.35) es, como de costumbre, el vector de fuerzas nodales
de equilibrio utilizado para el ensamblaje, con
(e)
qi = [P̄i , M̄xi , M̄yi ]T (8.38)
donde P̄i y M̄xi , M̄yi son la carga puntual vertical y los dos momentos flectores
que equilibran el nodo i.
La integración de los términos de la matriz de rigidez no es complicada y en
la Figura 8.7 se presenta su forma explı́cita para un elemento MZC homogéneo e
isótropo. En la Figura 8.8 se muestra la expresión del producto DBi para cálculo
de los momentos flectores. Finalmente, la expresión del vector de fuerzas nodales
equivalentes para una carga uniformemente repartida es
1 a b 1 a b 1 a b 1 a b
f (e) = 4qab , , , , − , , , − , − , , , − (8.39)
4 12 12 4 12 12 4 12 12 4 12 12
8.13
(e) (e) (e) (e) Et3
Figura 8.7
K(e) = D [K1 + K2 + K3 + K4 ] ; D=
12(1 − ν 2 )
6 6
8a2 0 0
6a 0 0
6b 0 8b2
0
−6 −6a 0 6 3 0 3b 6
6a 4a2 0 −6a 8a2 Simétri ca 0 0 0 0 0 Simétrica
0 0 0 0 0 0 3b 0 4b2 6b 0 8b2
(e) b (e) a
K1 = 3 ; K2 = 3
6a −3 −3a 0 3 −3a 0 6 6b −3 0 −3b −6 0 −6b 6
3a 2a2 0 −3a 4a2 0 −6a 8a2 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 3b 0 2b2 6b 0 4b2 −6b 0 8b2
3 3a 0 −3 3a 0 −6 6a 0 6
−6 0 −6b −3 0 −3b 3 0 −3b 6
3a 4a2 0 −3a 2a2 0 −6a 4a2 0 6a 8a2 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 6b 0 4b2 3b 0 2b2 −3b 0 4b2 −6b 0 8b2
8.14
1 21
0 8a2
ab 2ab 0
3a 0 8b2
3b
−1 0 −b 1 −21 −3a −3b 21
0 0 0 −a 0 3a −2a2 0 −3a 8a2
Simétrica Simétrica
−b 0 0 b −2ab 0 −3b 0 −8b2 3b 0 8b2
(e) ν ; K(e) = 1 − ν
K3 = 4
30ab
2ab 1 0 0 −1 a 0 1 21 3a 3b −21 3a −3b 21
0 0 0 a 0 0 −a 0 −3a 2a2 0 3a −8a2 0 −3a 8a2
0 0 0 0 0 0 −b 2ab 0 −3b 0 2b2 3b 0 −2b 2 −3b 0 8b2
−1 −a 0 1 0 0 −1 0 b 1
−21 −3a −3b 21 −3a 3b −21 3a 3b 21
Cálculo de Estructuras por el Método de Elementos Finitos
Mx
4
(e)
σ f = My
σ̂ = D̂f Bi ai
Mxy i=1
(d̂12 Nixx + dˆ12 Niyy ) dˆ11 N̄ixx ¯ yy
dˆ11 N̄ i
¯ yy
D̂f Bi = − (d̂21 Nixx + dˆ22 Niyy ) dˆ21 N̄ixx dˆ22 N̄ i
xy
2dˆ33 Nixy 2dˆ33 N̄ixy ˆ ¯
2d33 N̄ i
1 1
Nixx = − (3ξi ξ + 3ξi ηi ξη) N̄ixx = (3ξ + ξi ηi η + 3ηi ξη + ξi )
4a2 4a
1 1
Niy = − 2 (3ηi η + 3ξi ηi ξη) N̄ixy = (3ηi ξ 2 + 2ξi ηi ξ − ηi )
4b 8b
Nixy =
1 ¯ yy =
(4ξi ηi − 3ξi ηi ξ 2 − 3ξi ηi η 2 ) N̄
1
(3η + ξi ηi ξ + 3ξi ξη + ηi )
i
8ab 4b
¯ xy =
N̄
1
(3ξi η 2 + 2ξi ηi η − ξ)
i
8a
dˆij es el elemento ij de la matriz D̂f
Figura 8.8 Elemento de placa rectangular de cuatro nodos MZC. Matriz D̂Bi
para cálculo de momentos flectores.
8.15
Cálculo de Estructuras por el Método de Elementos Finitos
parcela obtenidos del cálculo, coinciden con los exactos correspondientes al campo
de movimientos prescrito [O3,Z3].
Por ejemplo, aunque es imposible garantizar la conformidad del elemento MZC
rectangular tomando como variables nodales la flecha y sus primeras derivadas,
el elemento satisface el criterio de la parcela, lo que asegura su convergencia al
disminuir el tamaño de la malla [Y1].
Desgraciadamente, para formas cuadriláteras arbitrarias deja de satisfacerse el
criterio de la parcela, perdiéndose todas las garantı́as de obtener convergencia.
Por consiguiente, el elemento cuadrangular de cuatro nodos MZC no es fiable para
usos prácticos. No obstante, en su forma rectangular es un elemento muy preciso
y puede utilizarse sin problemas [O3].
8.16
PLACAS DELGADAS Y GRUESAS
Los elementos triangulares son de gran interés práctico para análisis de placas
de formas irregulares. No obstante, en la formulación de dichos elementos se
encuentran las mismas dificultades para garantizar su conformidad que en los
elementos rectangulares. Presentaremos seguidamente los conceptos básicos de
algunos de los elementos de placa de Kirchhoff triangulares no conformes y
conformes más populares.
constante.
Tocher [T8] propone agrupar los términos a8 y a9 del polinomio cúbico como
8.17
Cálculo de Estructuras por el Método de Elementos Finitos
8.18
PLACAS DELGADAS Y GRUESAS
8.19
Cálculo de Estructuras por el Método de Elementos Finitos
siendo ȳ ortogonal al lado 23. Para los triángulos 412 y 431 se utilizan
aproximaciones similares. La omisión del término x̄2 ȳ garantiza que el giro normal
varı́a linealmente a lo largo de los contornos, mientras que la flecha varı́a de forma
cuadrática. Una vez obtenidas las matrices de rigidez de cada subelemento en ejes
locales, se efectúa el ensamblaje en ejes globales y se eliminan los grados de libertad
interiores, imponiendo al mismo tiempo continuidad de la derivada normal en los
puntos centrales de los lados interiores [C5,Y1].
8.20
PLACAS DELGADAS Y GRUESAS
El elemento anterior puede mejorarse añadiendo tres nodos laterales a los que
se asigna como variable la derivada normal ∂w ∂n [A7], [B4], [I3] (Figura 8.16b). En
dicho caso las funciones de forma son polinomios quı́nticos completos y ∂w ∂n varı́a
según un polinomio de cuarto grado a lo largo de los lados.
Más información sobre el desarrollo de estos elementos puede encontrarse en
[Y1] y [Z3,8] además de las referencias ya citadas. Pese a su excelente precisión los
elementos que incorporan las derivadas cruzadas como variables nodales no han
sido muy populares, debido fundamentalmente a las dificultades para su utilización
en análisis de láminas.
8.21
Cálculo de Estructuras por el Método de Elementos Finitos
8.22
PLACAS DELGADAS Y GRUESAS
8.23
Cálculo de Estructuras por el Método de Elementos Finitos
8.24
PLACAS DELGADAS Y GRUESAS
8.25
Cálculo de Estructuras por el Método de Elementos Finitos
de la placa, permanecen al deformarse sobre una misma recta, sin que ésta tenga
que ser necesariamente ortogonal a la deformada del plano medio (Figura 8.19).
El lector reconocerá la analogı́a de esta hipótesis con la establecida para el giro
de la sección en la teorı́a de vigas de Timoshenko. Esta analogı́a será de gran
utilidad para interpretar muchos aspectos comunes a ambas teorı́as.
donde θx y θy son los ángulos que definen el giro de la normal. Puede comprobarse
que el campo de desplazamientos anterior coincide con el expresado por la ecuación
(8.2) para la teorı́a de Kirchhoff. El vector de movimientos se define de igual forma
como
u = [w, θx, θy ]T (8.45)
8.26
PLACAS DELGADAS Y GRUESAS
∂u ∂θx
εx = = −z
∂x ∂x
∂v ∂θy
εy = = −z
∂y ∂y
∂w
εz = 0
∂z
∂u ∂v ∂θx ∂θy
γxy = + = −z +
∂y ∂x ∂y ∂x
∂u ∂w ∂w
γxz = + = −θx + = −φx
∂z ∂x ∂x
∂v ∂w ∂w
γyz = + = −θy + = −φy (8.48)
∂z ∂y ∂y
Se aprecia en (8.48) que la hipótesis de no ortogonalidad de la normal se traduce
en que las deformaciones transversales γxz y γyz no son nulas, siendo precisamente
su valor (absoluto) el de los giros φx y φy , respectivamente, que adquieren ası́ un
interesante significado fı́sico. Asimismo, se aprecia que dichas deformaciones (y
por consiguiente las respectivas tensiones) son independientes de la coordenada z.
Adviértase que la condición de deformaciones transversales nulas implica θx =
∂w y θ = ∂w , recuperándose la hipótesis de ortogonalidad de la normal de la
∂x y ∂y
teorı́a de Kirchhoff, como era de esperar.
8.27
Cálculo de Estructuras por el Método de Elementos Finitos
8.28
PLACAS DELGADAS Y GRUESAS
con G = 2(1+ν)E .
8.29
Cálculo de Estructuras por el Método de Elementos Finitos
σf
σ̂
+t zDf εf
2
σ =
σ̂ ······ = · · · · · · · · · · · · dz (8.55)
− 2t
σc
σ̂
Dcεc
y operando
−z ∂θx
∂x
∂θ
zDf −z ∂yy
σ
σ̂ f
+t
2
−z ∂θx + ∂θy
σ̂
σ = ······ = ∂y ∂x dz =
− 2t
σc
σ̂
························
∂w
− θ
∂x x
Dc ∂w
∂y − θ y
8.30
PLACAS DELGADAS Y GRUESAS
− ∂θ
∂x
x
t ∂θy
+ 2 2 − ∂y
z dz Df
3
− 2t
t
12 Df ε ε̂f
εf
D̂f ε̂
− ∂θx ∂θy
+ ······
······
= ∂y ∂x = =
······························ tDcε̂εc D̂cε̂εc
∂w
+2t
∂x − θ x
t dz Dc ∂w
−2 − θ
∂y y
(8.56)
t D3
donde D̂f = 12 f y D̂c = tDc (8.57)
..
σ̂ f
σ D̂ f . 0
ε̂
εf
σ =
σ̂ ······ =
· · · · · · · · · · · · · · · ······
= D̂ε̂ε (8.60)
.
σc
σ̂ 0 .. D̂c εε̂c
! "# $! "# $
D̂ ε
ε̂
8.31
Cálculo de Estructuras por el Método de Elementos Finitos
w
n Ni wi
u = θx = Ni θxi =
θy i=1 Ni θyi
8.32
PLACAS DELGADAS Y GRUESAS
w1
θ x
1
.. ..
θ y 1
0
N1 0 0 . . Nn 0 · · ·.· · ·
.. .. .
= 0 N1 0 . ······ . 0 Nn
0 . =
.. ..
· · · · · ·
0 0 N1 . . 0 0 Nn
wn
θx
n
θyn
a1
= [N1 , N2, . . . , Nn ] ... = Na(e) (8.63)
an
donde
(e)
a1
(e)
N = [N1 , N2 , . . . , Nn ] , a(e) = a2
..
.
(e)
an
y
Ni 0 0
wi
(e)
Ni = 0 Ni 0 , ai = θxi (8.64)
0 0 Ni θyi
8.33
Cálculo de Estructuras por el Método de Elementos Finitos
− ∂θ
∂x
x
− ∂N iθ
x i
∂x
∂θy
− ∂y
∂N
− i θ y i
ε
ε̂
∂y
f ∂θy n ∂Ni ∂Ni
∂θx
ε̂ε = · · · · · · = − ∂y + ∂x = − ∂y θ x i + ∂x i θ y =
··················
εε̂c
············
i=1
∂w
∂N i −
− θ x
∂x w i N i θ x i
∂x
∂w − θ ∂N iw − N θ
∂y y ∂y i i y i
(e)
n
B a1
(e)
= f i ai = [B1 , . . . , Bn ] .. = Ba(e)
Bci
.
i=1 (e)
an
(8.65)
0 − ∂N i 0
∂x ∂Ni
−Ni 0
− ∂N i ∂x
con Bfi = 0 0 ∂y y Bci = ∂Ni
0 −Ni
0 − ∂N
∂y
i − ∂N
∂x
i ∂y
(8.66)
(e)
BT D̂BdA a(e) − NqdA = q(e) (8.68)
A A(e)
8.34
PLACAS DELGADAS Y GRUESAS
o
K(e) a(e) − f (e) = q(e) (8.69)
donde
(e)
Kij = Bi T D̂Bj dA (8.70)
A(e)
(e)
fi = Ni [q, 0, 0]T dA (8.71)
A(e)
son la submatriz de rigidez que conecta los nodos ij y el vector de fuerzas nodales
equivalentes del nodo i debido a una carga repartida vertical. El vector de fuerzas
nodales de equilibrio del nodo i coincide con la ec.(8.38).
Haciendo uso de (8.60) y (8.66) se puede transformar la expresión de la matriz
de rigidez del elemento en la forma siguiente:
T BfT T
K(e) = (e)
[Bf , Bc ] D̂ dA =
A Bc
(e) (e)
= BTf D̂f Bf + BTc D̂c Bc dA = Kf + Kc (8.72)
A(e)
donde
(e)
Kf = BTf D̂f Bf dA (8.73a)
A(e)
y
(e)
Kc = BTc D̂c Bc dA (8.73b)
A(e)
8.35
Cálculo de Estructuras por el Método de Elementos Finitos
(Kf + Kc )a = f (8.77)
Et 3
Sacando factores comunes 12(1−ν2) y Gt en Kf y Kc , respectivamente (ver
ecs.(8.73)y (8.57)), se tiene
Et3
K̄ + GtK̄c a = f (8.78)
12(1 − ν 2 ) f
8.36
PLACAS DELGADAS Y GRUESAS
con
12(1 − ν 2 )G
α = (8.79b)
Et2
El segundo miembro de (8.79a) es, pues, del orden de magnitud de la solución
exacta para placas delgadas ak . Observando dicha ecuación vemos que para
t −→ 0, α −→ ∞. Por consiguiente, al hacerse la placa más delgada los términos de
cortante van progresivamente dominando la solución, de forma que la contribución
de Kf puede despreciarse. Ası́, cuando t −→ 0 la ec.(8.79a) tiende a
1
αK̄c a = O(ak ) y K̄c a = O(ak ) = 0 (8.80)
α
Se aprecia, pues, que en el lı́mite para α −→ ∞ se obtiene una solución
infinitamente más rı́gida que la correspondiente a la teorı́a de placas delgadas
y la única forma de obtener una solución diferente de la trivial a = 0 es que la
matriz Kc sea singular.
Similarmente al caso de flexión de vigas de Timoshenko, la integración reducida
de Kc proporciona la singularidad necesaria, y de nuevo hay que hacer uso de la
condición (7.68) para estudiar en cada malla si existe o no singularidad.
El proceso de subintegrar la matriz de rigidez de cortante manteniendo la
integración exacta de la matriz de rigidez de flexión se denomina integración
selectiva.
En la Figura 8.23 se muestran los elementos cuadriláteros y triangulares más
populares que presentan un comportamiento razonable para análisis de placas
gruesas y delgadas gracias a la integración selectiva.
Hay que señalar que en los elementos cuadriláteros de 4 y 9 nodos la integración
selectiva introduce dos y un mecanismos adicionales en la matriz de rigidez,
respectivamente. Estos mecanismos pueden propagarse en algunos casos, dando
lugar a soluciones incorrectas. Pese a ellos, su comportamiento es bueno en todo
el rango de espesores delgados y gruesos para la mayor parte de los problemas de
interés práctico [O3].
El cuadrilátero serendı́pito de 8 nodos con integración selectiva está libre
de mecanismos, aunque lamentablemente no da buenos resultados para placas
delgadas con relación de espesor/lado menores que 10−2 [O3].
El elemento triangular de 6 nodos con integración reducida para Kc presenta
un excelente comportamiento para análisis de placas delgadas y gruesas.
Desgraciadamente, el sencillo elemento triangular de tres nodos sufre de un alto
nivel de bloqueo, incluso con la integración reducidad de la matriz Kc .
En el apartado siguiente se describe una técnica alternativa, basada en la
selección “a priori” de un campo de deformaciones de cortante que puede satisfacer
las condiciones de Kirchhoff en el lı́mite de placas delgadas. Esta técnica de
deformaciones de cortante “impuestas” es muy popular hoy en dı́a y en base a
ella se ha desarrollado un elemento de placa cuadrilátero de 4 nodos que goza de
gran aceptación. En el Apartado 8.16 se presentan los detalles de la formulación
de este elemento.
8.37
Cálculo de Estructuras por el Método de Elementos Finitos
· · · · · · + αn (wi , θ i )ξ p ηq = 0 (8.81)
8.38
PLACAS DELGADAS Y GRUESAS
α j (wi , θ i ) = 0 ; j = 1, n (8.82)
La ec (8.82) impone una relación lineal entre los giros y flechas nodales que
generalmente puede también interpretarse desde un punto de vista fı́sico. Los
elementos que satisfacen (8.82) pueden reproducir de forma natural las condiciones
lı́mites de placa delgada sin que aparezca el efecto de bloqueo.
Sin embargo, en muchos elementos las αj son únicamente función de los giros
nodales y la condición αj (θθ i ) = 0 es demasiado restrictiva (e incluso no natural)
lo que provoca el bloqueo de la solución.
Ası́, por ejemplo, consideremos el sencillo elemento cuadrilátero de 4 nodos de
la Figura 8.23. De sus funciones de forma es fácil deducir que la deformación de
cortante γxz viene dada por (Capı́tulo 4)
4
∂w ξ 1 ξη η
γxz = − θx = ( i wi − θxi ) + ( i i wi − i θxi )η −
∂x i=1 4a 4 4a 4
ξ ξη (8.83)
− ( i θxi )ξ − ( i i θxi )ξη = α1 (wi , θxi ) +
4 4
+ α2 (wi , θxi )η + α3 (θxi )ξ + α4 (θxi )ξη
8.39
Cálculo de Estructuras por el Método de Elementos Finitos
(e)
donde γ k son los valores de las deformaciones de cortante transversal en m
puntos dentro del elemento y Nγ las correspondientes funciones de interpolación.
Combinando (8.85) y (8.65) puede escribirse
m
εε̂ = Nγk Bck a(e) = B̂c a(e) (8.86)
k=1
8.40
PLACAS DELGADAS Y GRUESAS
al propuesto por Bathe y Dvorkin [B1] y Hinton y Huang [H5] que se presenta en
el Apartado 8.16.
Es interesante advertir que el campo de deformaciones de cortante impuestas
para este caso conduce a la misma matriz de rigidez de cortante que la que se
obtendrı́a a partir del elemento original utilizando la cuadratura de integración
(e)
para Kc que se muestra en la Figura 8.24b.
Para un cuadrilátero general con coordenadas isoparamétricas ξ y η la situación
es idéntica en cuanto al primer punto anterior se refiere, si se utilizan γξ y γη en
lugar de γxz y γyz (por sencillez denominaremos γξ y γη a las deformaciones γξζ
y γηζ en el sistema de coordenadas normalizado).
Zienkiewicz et al. [Z4,5] han demostrado que el campo de deformaciones
de cortante impuestas debe satisfacer ciertas condiciones. Ası́, a partir de una
formulación mixta en la que se interpolan de forma independiente la flecha, los
giros y las deformaciones de cortante, las condiciones que deben satisfacer dichas
interpolaciones son:
nθ + nw ≥ nγ
(8.90)
nγ ≥ nw
8.41
Cálculo de Estructuras por el Método de Elementos Finitos
Etapa 1
Consideremos el elemento en el sistema de coordenadas natural ξ, η.
La etapa inicial es la interpolación de las deformaciones de cortante transversal
naturales en el sistema ξ, η. Ası́
α1
α2
· · · ξ p ηq | ···
γ
γξ 1 ξ η ξη 0 0 0 0
= γη = 0 0 0 0 ··· 0 | 1 ξ η ··· r s
ξ η .
. = Aα
.
αnγ
(8.91)
[Nota: Para simplificar la notación denominaremos γ al vector de deformaciones
de cortante transversal en vez de εε̂c . Asimismo, como εc = ε̂εc no procede distinguir
las deformaciones de cortante generalizadas].
Las deformaciones de cortante en el sistema cartesiano x y z se obtienen
directamente de (8.91) como
γxz
γ= = J−1γ (8.92)
γyz
donde βi es el ángulo que la dirección ξ¯i forma con la dirección natural ξ. El sentido
de las direcciones ξ¯i sobre los lados del elemento puede escogerse de acuerdo con
la numeración creciente de los nodos de cada lado en cuestión.
8.42
PLACAS DELGADAS Y GRUESAS
Etapa 2
α = γ ξ̄
P(ξi , ηi , βi )α (8.94)
n
donde γ ξ̄ = [γξ1¯, γξ2¯, · · · , γξ¯ γ ]T contiene los valores de la deformación de cortante
impuesta en los nγ puntos. De (8.94) puede obtenerse
α = P−1γ ξ¯ (8.95)
Etapa 3
γ
γ ξ¯ = T(βi )γ̂ (8.96)
n n
γ = [γξ1 , γη1 , γξ2 , γη2 · · · , γξ γ , γη γ ]T .
donde γ̂
γ = A P−1 T γ
γ̂ (8.97)
Etapa 4
8.43
Cálculo de Estructuras por el Método de Elementos Finitos
Etapa 5
Etapa 6
Este elemento fue desarrollado por Bathe y Dvorkin [B1], [D4] y Hinton y Huang
[H6] como particularización de elementos similares basados en la introducción de
modos de cortante auxiliares propuestos por Mac Neal [M1] y Hughes et al. [H8].
Posteriormente, Donea y Lamain [D2] presentaron una interesante reformulación
del elemento basándose en conceptos similares a los estudiados en el Apartado
8.15.1. El punto de partida en todos los casos es el clásico elemento cuadrilátero
de 4 nodos de clase Co . Ası́, la geometrı́a y el campo de movimientos se interpolan
con las funciones bilineales (ver Capı́tulo 4).
Siguiendo los razonamientos del Apartado 8.15.1 se deduce la siguiente
aproximación de las deformaciones de cortante transversal en el sistema de
coordenadas natural ξ, η (Figura 8.24a)
γξ = α1 + α2 η
(8.103)
1 η 0 0
γη = α3 + α4 ξ ; i.e. A =
0 0 1 ξ
8.44
PLACAS DELGADAS Y GRUESAS
Las αi se obtienen colocando las deformaciones naturales γξ¯ en los cuatro puntos
que se muestran en la Figura 8.25, siendo
1
γξ¯
1 −1 0 0 α1
2
1 0 1 0
0 1 γ¯ 1 −1 0
0 1 α2 0 1
= ξ
y P−1 = (8.105)
1 1 0 0 α3
γ 3¯
2 0 1 0 1
ξ
0 0 1 −1 α4
γ4 0 1 0 −1
! "# $ ξ ¯
P
8.45
Cálculo de Estructuras por el Método de Elementos Finitos
8.46
PLACAS DELGADAS Y GRUESAS
Como vimos en el Apartado 8.9, una de las técnicas para desarrollar elementos
de placa de Kirchhoff conformes es introducir ciertas modificaciones en los
elementos de placa de Reissner-Mindlin de manera que se satisfagan de forma
discreta sobre el elemento las condiciones de Kirchhoff. De ahı́ el nombre de
elementos DK (Discretos de Kirchhoff).
La idea de los elementos de placa DK es original de Wempner et al. [W3,4],
quienes la utilizaron como un método de evitar los requisitos de continuidad
C1 de los elementos de placa de Kirchhoff. Las condiciones de placa delgada
(γxz = γyz = 0) se imponen en puntos discretos de elementos de clase Co
formulados en base a la teorı́a de placas de Reissner-Mindlin. Pese a su cierta
antigüedad, este método sólo se popularizó como resultado del éxito del elemento
de lámina “semi-loof” de Irons [I6,7], considerado por algunos como un elemento
DK, y también gracias a una reinterpretación del método DK en base a su analogı́a
con el de deformaciones de cortante impuestas. Ası́, mientras en esta última técnica
se impone una determinada variación de las deformaciones de cortante sobre el
elemento, en los elementos DK dicha variación conduce a que dichas deformaciones
sean efectivamente nulas sobre el elemento.
En los últimos años se han propuesto varios elementos de placa DK.
Presentaremos aquı́ el elemento triangular de tres nodos por su gran versatilidad
para análisis de placas y láminas delgadas.
8.47
Cálculo de Estructuras por el Método de Elementos Finitos
donde Ni son las clásicas funciones de forma del elemento triangular de seis
nodos de clase Co (Capı́tulo 4).
3) Se impone una variación lineal del giro normal θn (ver Figura 8.9a) a lo largo
de los lados (tres condiciones), es decir
1
θnk = (θni + θnj ) (8.110)
2
donde k denota el nodo intermedio del lado ij.
8.18 CONCLUSIONES
8.48
PLACAS DELGADAS Y GRUESAS
8.49
CAPÍTULO 10
10.1 INTRODUCCIÓN
10.1
Cálculo de Estructuras por el Método de Elementos Finitos
10.2
LÁMINAS DE REVOLUCIÓN Y ARCOS
10.3
Cálculo de Estructuras por el Método de Elementos Finitos
siendo z la distancia entre los puntos O y P . Ası́ pues, de (10.1) y (10.2) se tiene
De (10.3)–(10.5) se deduce
u = uo − z θ
(10.6)
w = wo
Las ecs.(10.6) son la versión unidimensional del campo de desplazamientos en
láminas planas (ec.(10.1)). Es decir, el desplazamiento en la dirección tangente u
es suma del desplazamiento de membrana uo y del correspondiente desplazamiento
por el giro de flexión z θ. Por otra parte, el desplazamiento en dirección normal
es constante a lo largo del espesor.
De (10.6) se define el vector de movimientos locales de un punto como
uo
u = wo (10.7)
θ
Igualmente puede definirse un vector de movimientos globales u, que se
relaciona con u por la transformación
uo cosφ senφ 0 uo
u =
w
o = −senφ cosφ 0 wo = L u (10.8)
θ 0 0 1 θ
De (10.3) y (10.4) se desprende que los movimientos uo , wo y θ son función
únicamente de la longitud de arco s. Esto es de gran importancia para la obtención
de las deformaciones [O3].
10.4
LÁMINAS DE REVOLUCIÓN Y ARCOS
10.5
Cálculo de Estructuras por el Método de Elementos Finitos
donde
E 1 ν αE
Dm = y Dc = (10.15)
1−ν ν
2 1 2(1 + ν)
10.6
LÁMINAS DE REVOLUCIÓN Y ARCOS
Nx σx
σ
σ̂ m
Nϕ σϕ
σm
...
...
...
t
t
...
2 2
σ =
σ̂
σf
σ̂
= Mx = t
z σ
x dz = t
z σ
m dz
(10.16)
−2 −2
...
Mϕ
z σϕ
...
σc
σ̂
. . .
. . .
σ
c
Qz τx z
σ m = [Nx , Nϕ ]T , σ̂
siendo σ̂ σ f = [Mx , Mϕ ]T y σ̂c = Qz .
Utilizando (10.14) y (10.16) y siguiendo un proceso similar al del Apartado
9.3.3, se encuentra la relación entre esfuerzos y deformaciones generalizadas locales
σ = D̂ εε̂
σ̂ (10.17)
siendo
ε
ε̂m
ε̂ε =
εε̂f
(10.18)
ε̂c
σ m = D̂m ε̂εm
σ̂ , σ f = D̂f ε̂εf
σ̂ y σ̂c = D̂c ε̂c (10.20)
t3
D̂m = tDm , D̂f = D , D̂c = tDc (10.21)
12 m
10.7
Cálculo de Estructuras por el Método de Elementos Finitos
δεεT σ dV
ε σ = δuT b t dA + δuT t dA + δuT
i pi dΓ (10.22)
V A A i Γ
dA = dx dz ds dz (10.24)
t
δεεT σ x dA = 2π δε̂εT εT
2
2π m σ m + δε̂ f z σ m + δ ε̂c σc x ds dz =
A s −t
2
= 2π δε̂εT
m σ σ̂ m + δε̂εT
f σ̂ f + δ ε̂c σ̂c x ds =
σ
s
= 2π δε̂εT σ̂ σ x ds (10.25)
s
10.8
LÁMINAS DE REVOLUCIÓN Y ARCOS
10.9
Cálculo de Estructuras por el Método de Elementos Finitos
10.10
LÁMINAS DE REVOLUCIÓN Y ARCOS
n
(e)
εε̂ = Bi ai = B a
(e)
(10.30)
i=1
siendo
y
∂Ni
∂s 0 0
N cosφ N senφ
i − ix 0
x
Bmi
—
—————— — — — — —
− − −
Bi =
Bfi = 0 0 −
∂Ni (10.32)
∂s
− − −
Ni cosφ
Bci 0 0 − x
— —————— — — — — —
∂Ni
0 ∂s −Ni
donde
(e)
Kij = 2π BT
i D̂ Bj x ds (10.34)
l(e)
es una submatriz tı́pica de la matriz de rigidez del elemento en ejes locales; f (e)
es el vector de fuerzas nodales equivalentes locales del elemento, cuya expresión en
ejes globales se dará más tarde; q(e) es, como de costumbre, el vector de fuerzas
nodales de equilibrio y l(e) la longitud del elemento.
10.11
Cálculo de Estructuras por el Método de Elementos Finitos
con
(e) (e)
Kij = LT Kij L (10.37)
donde
L 0 1
L 2
T(e) = ... .. (10.38)
.
0 L n
10.12
LÁMINAS DE REVOLUCIÓN Y ARCOS
Bmi
∂Ni ∂Ni
0 ∂s cos φ ∂s senφ
Ni
0 0 x
B
mi — — — — — — — — — — —
Bi = Bfi = Bi L =
0 0 − ∂Ni
∂s
Bfi
Bci
0 0 − Ni cos φ
x
— — — — — — — — ———
− ∂N i
∂s senφ
∂Ni
∂s cos φ −Ni
Bci
(10.39)
De esta forma la matriz de rigidez global se obtiene por
(e) (e) (e) (e) (e)
Kij = 2π BTi D̂ Bj x ds = Kmij + Kf + Kcij + Kmf (10.40)
l(e) ij ij
donde las matrices de rigidez anteriores se obtienen utilizando las nuevas Bmi , Bfi
y Bci de (10.39) en las expresiones (10.35).
La matriz de deformación Bi permite calcular directamente los esfuerzos locales
a partir de los movimientos globales por (9.45).
La ecuación de equilibrio en ejes globales se escribe ahora
10.13
Cálculo de Estructuras por el Método de Elementos Finitos
10.14
LÁMINAS DE REVOLUCIÓN Y ARCOS
+1
nq
D̂ D̂ Bj xJ (e) Wq =
(e)
Kij = 2π BTi Bj xJ (e) dξ = BTi
−1 q=1
q
nm ! " nf
(e) (e)
= Im Wqm + If Wqf +
qm
qf
qm =1 qf =1
! " n mf
(e)
nc
(e)
+ Ic Wqc + Imf W (10.45)
qc qmf qmf
qc =1 qmf =1
donde
nq nq
(e) (e)
T
fi = 2π Ni b(e) xtJ (e) Wq + 2π NTi t(e) xJ (e) Wq + 2πxi pi (10.47)
q q
q=1 q=1
10.15
Cálculo de Estructuras por el Método de Elementos Finitos
(e) ¯
Kij = 2π B̄Ti D̂ B̄j x̄ l(e) (10.48)
donde (¯·) indica valores en el centro del elemento (ξ = 0). La matriz B̄i se obtiene
(−1)i
directamente haciendo ∂N i
∂s = l(e) , Ni = 1/2 y x = x̄ en la ec.(10.39). La forma
(e)
explı́cita de Kij para este elemento para el caso de acoplamiento membrana-
flexión nula (Dmf = 0) puede verse en la Figura 10.8.
Por ser de gran interés práctico, presentamos seguidamente la expresión exacta
del vector de fuerzas equivalentes nodales del elemento troncocónico de dos nodos
para el caso de cargas másicas y repartidas constantes sobre el elemento. Esta es
10.16
Figura 10.8
(−1)i+j 2 (−1)i+j
(C d11 + S 2 d55 ) + d222 + | SC(d11 − d55 )+ |
(l(e) )2 4x (l(e) )2
i
(−1)j (−1)
| + S d21 | S d
+ C(e) ((−1)i d12 + (−1)j d21 ) (e) 55
2x̄l 2x̄l(e) 2l
| |
10.17
— — — — — — — — — — — — — — — — — — — — — — — ——————————————
| |
i+j i i
(e) (e) (−1) (−1) (−1)i+j 2 (−1)
Kij = C (l(e) )2 SC(d11 − d55 ) + 2x̄l(e) S d12 | (S d11 + C 2 d55 ) | − (e) C d55
i,j=1,2 (l(e) )2 2l
| |
— — — — — — — — — — — — — — — — — — — — — — — ——————————————
LÁMINAS DE REVOLUCIÓN Y ARCOS
| |
10.18
LÁMINAS DE REVOLUCIÓN Y ARCOS
wi
ai =
θi
(−1)i+j (−1)i
(l(e) )2 d33 − d
2l(e) 33
i+j
2πl x̄
(e)
+ d22
Kij = (e) (−1) i, j = 1, 2
(−1)j d11 +
−
(l(e))2 4x¯2
d
2l(e) 33
+ d12(e) [(−1)i + (−1)j ] + d433
2x̄l
D̂f 0
dij : componentes de D̂ = (ver (10.21)) evaluada en el centro del elemento
0 Dc
Figura 10.9 Matriz de rigidez del elemento de placa de revolución de dos nodos
de Reissner-Mindlin con integración reducida de un solo punto.
10.19
Cálculo de Estructuras por el Método de Elementos Finitos
10.20
LÁMINAS DE REVOLUCIÓN Y ARCOS
rectos, R = ∞, y
∂uo
alargamiento
∂s
ε̂ε = − ∂θ curvatura (10.58)
∂s
deformaciones de
∂wo
flexión
−θ cizalladura
∂s
Puede apreciarse que εε̂ se compone ahora de las deformaciones de flexión
del elemento de viga de Timoshenko del Capı́tulo 7 y la deformación axil
(alargamiento) del elemento de barra a tracción del Capı́tulo 2, que actúan de
forma desacoplada a nivel local . El acoplamiento entre los efectos axiles y de
flexión se produce en ejes globales al ensamblar las contribuciones de los diferentes
elementos.
Discretizando el arco en elementos rectos de n nodos se obtiene la matriz de
deformación generalizada local, como
n
(e)
ε̂ε
ε = Bi ai = B a(e) (10.59)
i=1
con
∂Ni
∂s
0 0
uoi
(e)
Bi =
0 0
−∂Ni
; ai = w (10.60)
∂s
oi
∂Ni θi
0 ∂s −Ni
10.21
Cálculo de Estructuras por el Método de Elementos Finitos
(e) (e)
donde KBT y KV T son las matrices de rigidez de los elementos de barra
ij ij
a tracción y de viga de Timoshenko estudiados en los Capı́tulos 2 y 7,
respectivamente.
La transformación de la matriz de rigidez local a ejes globales para el ensamblaje
sigue exactamente las mismas etapas descritas para láminas de revolución,
obteniéndose,
Kij = (e) BTi D̂ Bj dx
(e)
(10.64)
l
(e) ¯
Kij = B̄Ti D̂ B̄j l(e) (10.65)
donde (·) indica de nuevo valores en el centro del elemento. En la Figura 10.11 se
(e)
muestra la expresión de Kij de (10.65).
Figura 10.11 Matriz de rigidez del elemento de arco plano de dos nodos con
integración reducida uniforme de un punto.
10.22
LÁMINAS DE REVOLUCIÓN Y ARCOS
10.23
CAPÍTULO 9
ANÁLISIS DE LÁMINAS
CON ELEMENTOS PLANOS
9.1 INTRODUCCIÓN
Las estructuras laminares son muy comunes en numerosos campos de la
ingenierı́a. Como ejemplos podrı́amos citar: puentes, cubiertas, depósitos, cascos
de barco, fuselaje de aviones, carrocerı́as de vehı́culos, etc.
Tipológicamente las láminas pueden considerarse una generalización de las
placas al caso de superficie media no plana. Es precisamente esta no coplanaridad
la que confiere el carácter resistente de las láminas al permitir la aparición de
esfuerzos axiles (esfuerzos de membrana) que, juntamente con los de flexión,
contribuyen a dotar a las láminas de una capacidad portante muy superior a la de
las placas.
En general, podemos decir que las láminas son a las placas, lo que los arcos (o
las estructuras reticulares) son a las vigas. Por tanto, un buen conocimiento de
la influencia del axil en arcos y pórticos favorecerá sin duda la comprensión del
funcionamiento resistente de las estructuras laminares. En la Figura 9.1 se muestra
un sencillo esquema de la contribución de los esfuerzos axiles a la resistencia de
un pórtico y de una lámina formada por ensamblaje de dos placas.
Desde el punto de vista geométrico las láminas se clasifican por la forma de su
superficie media [T5], [V1]. Estudiaremos aquı́ los casos de láminas de superficie
media arbitraria (Capı́tulo 9) y las láminas de revolución (Capı́tulo 10).
La obtención de las ecuaciones de una lámina (equilibrio, constitutivas y
cinemáticas) es complicada, debido precisamente a la curvatura de su superficie
media [F3], [K2], [N3], [T5], [V1]. Una de las maneras más sencilas de sortear
este problema es estudiar el comportamiento de una lámina como si estuviese
compuesta de elementos planos de tamaño pequeño.
Parece intuitivo que la aproximación de la geometrı́a real será tanto más exacta
cuanto más pequeño sea el tamaño de la discretización utilizada, análogamente
al proceso de aproximación de una lı́nea curva por rectas progresivamente más
pequeñas (Figura 9.2). La idea anterior es la base de la teorı́a de elementos finitos
de “lámina plana” que se trata en este capı́tulo. Esta teorı́a es de gran interés
no sólo para estudiar de forma aproximada láminas de superficie media curva,
sino también como método natural de análisis de numerosas estructuras laminares
compuestas por elementos de placa ensamblados en el espacio. Ejemplos de estas
9.1
Cálculo de Estructuras por el Método de Elementos Finitos
9.2
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
9.3
Cálculo de Estructuras por el Método de Elementos Finitos
9.4
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
∂u −z ∂xx
∂θ
∂uo
∂x
∂x
εx
∂v
∂θ
∂vo
−z ∂y y
εy
∂y
∂y
∂u ∂v
∂θ
γ ∂θx y
∂y + ∂x
ε = x y = =
∂uo
+
∂vo
+ −z ( ∂y + ∂x )
······
∂y ∂x
······
·········
·········
γx z
∂u + ∂w
∂w o − θ
γy z
∂z ∂x
0
x
∂x
∂v + ∂w
0
∂w
o − θ
∂z ∂y y
∂y
(9.3)
o bien
ε̂εm
ε̂
zε f
ε = ······ + ······ (9.4)
0
ε̂εc
9.5
Cálculo de Estructuras por el Método de Elementos Finitos
donde
T
∂uo ∂vo ∂uo ∂vo
ε̂εm = , , ( + ) (9.5a)
∂x ∂y ∂y ∂x
T
∂θ ∂θy ∂θ ∂θy
εε̂f = − x , − , −( x + ) (9.5b)
∂x ∂y ∂y ∂x
T
∂w ∂w
ε̂εc = ( o − θx ), ( o − θy ) (9.5c)
∂x ∂y
son, respectivamente, los vectores de deformaciones generalizadas de membrana
(alargamientos), flexión (curvaturas) y cortante (cizalladuras). Adviértase que,
similarmente al caso de placas, las deformaciones γx z y γy z representan los giros
φx y φy , respectivamente, tal y como se puede deducir de (9.3) y la Figura 9.5.
Por otra parte, de las ecs. (9.3)-(9.5) se deduce que:
a) la deformacion total de un punto se obtiene sumando las deformaciones de
membrana a las del estado de placa, lo cual es consecuencia directa del campo
de desplazamientos escogido.
b) los vectores de deformaciones generalizadas de membrana y de cortante
contienen las deformaciones “en el plano” y transversales al mismo,
respectivamente, y
c) el vector de deformaciones generalizadas de flexión contiene las tres curvaturas
del plano medio.
σx
εx
σy
..
εy
σ 0
τx y f Df . γ
σ =
= ··· = ··· ··· xy
··· ······ = Dε (9.6)
······
..
τ
σc 0 . Dc
γ
x z
x z
τy z γy z
donde para material isótropo
1 ν 0
E
Df = ν 1 0
1 − ν2 1−ν
0 0 2
1 0 E
Dc = αG con G= (9.7)
0 1 2(1 + ν)
9.6
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
9.3.4 Esfuerzos
El vector de esfuerzos locales en un punto del plano medio se define por
Nx
σx
Ny
σy
Nx y
τ
σ
σ̂
x y
σf
···
m
···
···
Mx
+t
z σx
+t ···
σ = σ̂ dz
2 2
σ̂ σ f = = t z σ
dz = z σ f (9.10)
M
− − t
···
y 2
y 2 ···
M y
z τ
σ
x
x y
σc
σ̂ c
···
· · ·
τx z
Qx
Qy τy z
σ m , σ
donde σ̂ σ̂ f y σ̂
σ c son los vectores de esfuerzos locales de membrana, flexión
y cortante, respectivamente, y t el espesor. Los esfuerzos de flexión y cortante
coinciden con los estudiados para placas. El vector de esfuerzos de membrana lo
forman los tres axiles Nx , Ny y Nx y contenidos en el plano medio. El convenio
de signos se muestra en la Figura 9.6.
La relación entre esfuerzos y deformaciones generalizadas locales se obtiene
combinando (9.8)–(9.9) y (10.10) como
σm
σ̂
σ f
···
+t
···
σ = dz =
2 σ
σ̂ σf
σ̂ = z
− 2t
f
···
···
σ
σ̂ c σc
9.7
Cálculo de Estructuras por el Método de Elementos Finitos
Df (ε̂εm + z ε̂εf )
εε̂m
+t
············
···
2
=
z D (ε̂
ε
m + z ε
ε̂
f )
dz = D̂
ε
ε̂
f
= D̂ε̂ε (9.11)
− 2t f
············
···
Dcε̂εc ε̂εc
σ y el de
siendo D̂ la matriz constitutiva que relaciona el vector de esfuerzos σ̂
deformaciones generalizadas ε̂ε en ejes locales. De (9.11) se deduce
+t Df z Df 0 D̂m D̂mf 0
2
D̂ = t
z D
f z 2 Df 0
dz = D̂mf D̂f 0
(9.12)
−2
0 0 Dc 0 0 D̂c
con
+t +t
D̂m Df dz D̂mf
2
z Df dz
2
= ; =
− 2t −2t
+t +t (9.13)
2 2 2
D̂f = t
z Df dz ; D̂c = Dc dz
−2 − 2t
donde D̂m , D̂f y Dc son las matrices constitutivas generalizadas correspondientes
a esfuerzos de membrana, flexión y cortante, respectivamente, y D̂mf es la matriz
constitutiva de acoplamiento membrana-flexión.
Las expresiones anteriores son válidas para el caso más general en el que
las propiedades del material estén hetereogéneamente distribuidas a través del
espesor (por ejemplo, en el caso de láminas formadas por materiales compuestos
de propiedades mecánicas diferentes, o bien en láminas de hormigón armado con
distribuciones no simétricas de armaduras). Es fácil advertir que si existe simetrı́a
de las propiedades del material con respecto al plano medio, o bien el material
es homogéneo, D̂mf = 0 y cada vector de esfuerzos se puede calcular de forma
desacoplada a partir de sus correspondientes deformaciones generalizadas por
σ̂ m = D̂mε̂εm
σ ; σ f = D̂f ε̂εf
σ̂ ; σ̂ c = D̂cε̂εc
σ (9.14)
9.8
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
9.9
Cálculo de Estructuras por el Método de Elementos Finitos
9.10
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
donde B y Bi son las matrices de deformaciones generalizadas locales del elemento
y de un nodo i, respectivamente. Esta última puede escribirse como
Bmi
Bi =
Bf (9.22)
i
Bci
9.11
Cálculo de Estructuras por el Método de Elementos Finitos
0 0 0 − ∂N
∂x
i 0
i
Bfi = 0 0 0 0 − ∂N
∂y
(9.24)
∂Ni
0 0 0 − ∂y − ∂Ni
∂x
∂Ni
0 0 ∂x −Ni 0
Bci = ∂Ni
(9.25)
0 0 ∂y 0 −Ni
donde T
t = tx , ty , tz , mx , my (9.27)
siendo Rx , Ry y Rz las fuerzas puntuales que actúan en el nodo i del elemento
i i i
según direcciones x , y , z , respectivamente, y Mx , Mz los momentos nodales
i i
contenidos en los planos x z e y z .
Operando en la forma usual puede obtenerse la ecuación matricial de equilibrio
de un elemento aislado por
9.12
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
(e)
fi = NTi ti dx dy (9.30b)
A(e)
donde
(e)
B mi D̂m Bmj dx dy
T
Kmij =
A(e)
(e)
Kf = BfTi D̂f Bfj dx dy
ij
A(e)
(e)
Kcij = BcTi D̂c Bcj dx dy
A(e)
(e) (e) T
Kmf = (e)
Bm
T
i
D̂mf Bfj dx dy = Kf m (9.32)
ij A ij
9.13
Cálculo de Estructuras por el Método de Elementos Finitos
(e) (e)
donde KT P y KRM son las matrices de rigidez de los elementos de tensión plana
(Capı́tulo 4) y de placa de Reissner-Mindlin (Capı́tulo 9) con la misma tipologı́a
y número de nodos que el elemento de lámina plana utilizado.
Por consiguiente, si no existe acoplamiento membrana-flexión, la matriz de
rigidez local de un elemento de lámina plana puede obtenerse directamente
ampliando la matriz de rigidez para el caso de flexión de placas con la del elemento
de tensión plana correspondiente. En terminologı́a de cálculo de estructuras
podrı́amos decir que a nivel local los esfuerzos de membrana equilibran las acciones
contenidas en el plano del elemento, mientras que las acciones normales a dicho
plano provocan un estado de flexión independiente, pudiendo obtenerse, siempre
a nivel local, los movimientos, deformaciones y tensiones de ambos estados de
manera totalmente desacoplada. El “acoplamiento” entre los estados de membrana
y flexión se produce al ensamblar en ejes globales las ecuaciones de rigidez, tema
que se trata en el apartado siguiente.
9.14
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
donde
(e) T
ai = uoi , voi , woi , θxi , θyi , θzi
(e) T
fi = fxi , fyi , fzi , Mxi , Myi , Mzi (9.35)
donde
λx x λx y λx z (e)
λ (e) = λy x λy y λy z
(9.37)
λz x λz y λz z
9.15
Cálculo de Estructuras por el Método de Elementos Finitos
De (9.34) se deduce
donde
1 2···n
1
L(e)
T(e) = ... 2 (9.40)
..
5n × 6n .
L(e)
n
T T
q(e) = T(e) q(e) = T(e) K(e) a(e) − f (e) =
T T
= T(e) K(e) T(e) a(e) − T(e) f̄ (e) = K(e) a(e) − f (e) (9.41)
T
Bi D̂ Bj L(e) dx dy
(e) (e) T
Kij = L (e)
=
A
= BTi D̂ Bj dx dy (9.43)
A(e)
donde
Bi = Bi L(e) (9.44)
9.16
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
n
n
(e)
σ = D̂εε̂ = D̂ Bi ai = D̂ Bi L(e) ai
(e)
σ̂ =
i=1 i=1
n
= D̂ Bi ai = D̂ Ba(e) (9.45)
i=1
Tanto la matriz de rigidez global K(e) de (9.43) como el vector de fuerzas global
f (e) se calculan haciendo uso de integración numérica. Para ello, es esencial definir
primeramente las coordenadas de los nodos del elemento con respecto a los ejes
locales x y z . Esto puede hacerse mediante una transformación de coordenadas
idéntica a la (9.34) para los desplazamientos. Ası́, suponiendo que los orı́genes de
los sistemas local y global coinciden, se tiene
9.17
Cálculo de Estructuras por el Método de Elementos Finitos
siendo
9.18
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
y el vector unitario es
(e)
λx x
x
(e) 1 ij
vx = λx y = y
(e) ij
(9.51b)
λx z lij zij
siendo
(e)
lij = (x2ij + yij
2 + z 2 )(e)
ij (9.52)
(e)
λz x y z − zij yim (e)
(e) 1 (e) (e) 1 ij im
vz = λ = (Vij ∧ Vim ) = (e) xim zij − zim xij
zy (e) (e)
|Vij ∧ Vim | dz xij yim − yij xim
λz z
(9.53)
9.19
Cálculo de Estructuras por el Método de Elementos Finitos
(e) (e)
dz = (yij zim − zij yim )2 + (xim zij − zim xij )2 + (xij yim − yij xim )2
(e)
Es fácil deducir que en un elemento triangular, dz representa el doble de su
área, lo que puede simplificar los cálculos.
Finalmente, los cosenos directores del eje y se obtienen como producto vectorial
de los vectores unitarios en direcciones z y x . Ası́
λ (e) λ λ − λx y λz z (e)
y x
zy xz
(e) (e) (e)
vy = λy y = vz ∧ vx = λx x λz z − λz x λx z (9.54)
λy z λz x λx y − λz y λx x
9.20
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
9.21
Cálculo de Estructuras por el Método de Elementos Finitos
9.22
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
9.23
Cálculo de Estructuras por el Método de Elementos Finitos
∂w ∂w
u = uo − z , v = vo − z , w = wo (9.57)
∂x ∂y
donde las matrices D̂m , D̂mf y D̂f se obtienen por (9.13). De nuevo, para material
homogéneo, o con propiedades mecánicas simétricas con respecto al plano medio,
9.24
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
donde
T
∂wo ∂wo
u = uo , vo , wo , ,
∂x ∂y
T
(e) ∂wo ∂wo
ai = uoi , voi , woi , ( )i , ( )i (9.64)
∂x ∂y
9.25
Cálculo de Estructuras por el Método de Elementos Finitos
..
Ni 0 . 0 0 0 m ..
.. Ni . 0
0 Ni . 0 0
0 =
..
Ni =
··· .. ······ . ······
(9.65)
··· . ··· ··· ··· .. f
.. 0 . Ni
0 0 . Pi P̄i P̄¯i
9.26
ANÁLISIS DE LÁMINAS CON ELEMENTOS PLANOS
9.27
Cálculo de Estructuras por el Método de Elementos Finitos
9.12 CONCLUSIONES
9.28