Capítulo 12: El Método de Elementos Finitospolisoft.cl/clubpenrose/data/tg/tg_metodoelementosfinitos.pdf · El método de elementos finitos es una técnica de discretización física

Análisis Estático y Dinámico de Estructuras Tomás Guendelman Capitulo 12: El Método de Elemento Finitos

1

Capítulo 12: El Método de Elementos Finitos

INTRODUCCION

El método de elementos finitos es una técnica de discretización física de sólidos deformables que permite determinar en forma aproximada el estado de tensiones y deformaciones que se produce por la acción de solicitaciones externas. La discretización física consiste en el reemplazo del sólido original por un modelo físico equivalente, regido por ecuaciones de comportamiento que poseen ventajas relativas frente al sólido original (por ejemplo, ecuaciones diferenciales se transforman en ecuaciones algebraicas). En forma alternativa a la discretización física existe la discretización matemática, que consiste en la solución aproximada de las ecuaciones del problema original, mediante el empleo de técnicas numérica (por ejemplo, métodos de diferencias finitas, de mínimos cuadrados, Galerkin, etc.). El método de elementos finitos consiste en la división del sólido en un número finito de pequeños trozos o elementos que se conectan en sus nudos. Para estos pequeños elementos se introducen simplificaciones en sus relaciones cinemáticas de modo que los desplazamientos en cualquiera de sus puntos puedan ser expresados en función única de los desplazamientos de sus nudos. Se ensamblan posteriormente todos los elementos finitos entre sí mediante la aplicación de las condiciones de compatibilidad cinemática y de equilibrio en los nudos, generándose un sistema discreto (número finito de valores).


2

Consideremos un sólido plano que dividimos en elementos finitos triangulares, como muestra la figura.

x

y

0

y,v

x,u

a c

b

va

ua

vc

v(x,y) uc

u(x,y

vb

ub

0

( ){ } ( )

( )

=× y,xv

y,xu)y,x(r

12

( ){ }6 1×

=

q

uvuvuv

a

a

b

b

c

c

( ){ }

( )( )[ ]

( ){ }( )

r x y x y q, ,=× × ×2 1 2 6 6 1

θ (1)

La discretización del sólido se muestra en la primera figura, en la que podemos apreciar que todos los elementos son triangulares. De esta discretización extraemos un elemento cualquiera que se observa ampliado en la segunda figura.


3

Por tratarse de un problema de elasticidad plana, el movimiento de este elemento finito queda descrito totalmente por dos funciones continuas u(x,y) y v(x,y). Si se introducen simplificaciones de carácter cinemático, como se ha indicado anteriormente, estas funciones continuas pasan a tener expresiones que quedan completamente definidas por un número discreto de parámetros denominados grados de libertad. Si se escoge como grados de libertad los seis desplazamientos nodales de este elemento finito, las aproximaciones que introduzcámos pueden tomar la forma de la expresión (1). En ella se observa que conocida la matriz [θ(x,y)], denominada matriz de interpolación, los desplazamientos para cualquier punto del triángulo quedan totalmente definidos en función de sus seis desplazamientos nodales. Debe notarse que la matriz de interpolación es válida para un elemento y que las coordenadas (x,y) quedan restringidas al dominio geométrico correspondiente. El método de elementos finitos se sustenta en este punto de partida, a la vez que exige que la relación (1), denominada relación de comportamiento, cumpla con dos requisitos escenciales:

a) La matriz de interpolación debe ser capaz de reproducir movimientos de cuerpo rígido. Esto significa que en el elemento finito no se introduzcan restricciones geométricas producto de condiciones de apoyo implícitas.

b) Pese a que los elementos finitos se conectan entre sí sólo en sus puntos

nodales, es importante que implícitamente se genere compatibilidad geométrica en los bordes comunes entre elementos contiguos.

d

c

b a

posición original

d′

c′

b′ a′

posición deformadaINCOMPATIBLE

d′

c′

b′ a′

posición deformadaCOMPATIBLE


4

Se puede demostrar que la convergencia que experimentan las soluciones de problemas de elasticidad, en función del refinamiento de la malla de elementos finitos, tiene las siguientes formas:

refinamiento(por ejemplo:Nº de E. F. de la malla)

Solución(por ejemplo desplazamientosen un punto)

E.F. incompatibles

soluciónexactaE.F. compatibles

Ejemplo: P

δ

Refinamiento Al observar las curvas de convergencia podemos apreciar dos hechos importantes:

− Los modelos pierden rigidez en la medida que aumenta el refinamiento de la malla.

− Los elementos finitos compatibles convergen monotónicamente a la solución exacta.

La convergencia monotónica es muy importante pues permite, mediante un análisis de sensibilidad, medir el grado de precisión alcanzado. Pese al análisis recién hecho en función de la convergencia, hay muchos elementos que se aceptan en la práctica profesional y que son incompatibles. La razón para ello es que el grado de oscilación de su curva de convergencia sea muy pequeño y no constituya un problema mayor, al mismo tiempo que sus funciones de comportamiento tengan ventajas analíticas de importancia. En resumen, la primera de las dos condiciones que hemos impuesto a las funciones de comportamiento (capacidad de movimiento de cuerpo rígido) es obligatoria, mientras que la segunda es altamente recomendable, pudiendo ser eventualmente ignorada.


5

ELEMENTO FINITO TRIANGULAR

y, v

x, u

a

b

(xc , yc)

(xb , yb)

(xa , ya)

vc

uc

vb

va

ua

ub

v(x,y)

P(x,y)

u(x,y)

0

c

FUNCIONES DE COMPORTAMIENTO

Supongamos: ( ) yCxBAy,xu 111 ++= ( ) yCxBAy,xv 222 ++= Para (x, y) en el dominio del elemento finito. Justificación de las funciones escogidas:

a). Las funciones escogidas deben quedar plenamente definidas si se conocen los desplazamientos de los puntos nodales del elemento, que es la hipótesis fundamental del método de elementos finitos, que establece que el elemento se vincula al resto de la estructura sólo en sus puntos nodales. Por este motivo las funciones u(x,y) y v(x,y) poseen tres parámetros indeterminados cada una, denominados coordenadas generalizadas, que mediante transformaciones de tipo geométrico se pueden substituir por los desplazamientos u y v de los tres puntos nodales.

b). Se ha escogido un polinomio completo de primer grado en lugar de otra

expansión de tres parámetros indeterminados diferente a ésta. Demostraremos más adelante que es una elección satisfactoria pues cumple con las dos condiciones que validan el método. Por ahora hay dos factores previos que justifican esta elección. El primero es de carácter intuitivo y el segundo está relacionado con una simetría de contribución de las coordenadas x e y, lo que hace válida la elección, cualquiera sea la orientación del elemento.


6

A los dos factores anteriores debemos agregar la comprobación del cumplimiento de las dos condiciones que originan un elemento finito compatible:

- Capacidad de movimientos como cuerpo rígido.

- Compatibilidad geométrica inter-elemental. Lo que conduce a tres componentes independientes para los vectores {τ} {ε} y a una matriz constitutiva de (3×3). Las ecuaciones del movimiento de un cuerpo rígido son:

v0

φ0

u0

v(x, y)

u(x, y)P(x, y)

y

x

( )u x y u y, = −0 0φ ( )v x y v x, = +0 0φ Comparando estas ecuaciones con las funciones de comportamiento escogidas, observamos que ambas coinciden si se cumple que:

A u A vB BC C

1 0 2 0

1 2 0

1 0 2

00

= == == − =

;;;

φφ

Lo que permite concluir que las funciones de comportamiento escogidas tienen la capacidad de representar el movimiento de cuerpo rígido. Con respecto a la compatibilidad geométrica inter-elemental, podemos observar que si las funciones de comportamiento son de primer grado, los bordes rectos (y cualquier recta interior) del elemento finito triangular permanecerán siendo rectos en la posición deformada. Si dos elementos contiguos tienen compatibilidad geométrica explícita en sus dos nudos comunes, necesariamente el borde común tendrá una posición final idéntica para ambos elementos.


7

MATRIZ DE INTERPOLACION ( )[ ]θ x y,

Escribiendo en forma matricial las ecuaciones de comportamiento se tiene:

( )( )

u x y

v x y

x yx y

ABCABC

,

,

=

1 0 0 00 0 0 1

1

1

1

2

2

2

Simbólicamente: ( ){ }

( ) ( )( )[ ]{ }

( )r x y F x y, ,

2 1 2 6 6 1×

=

× ×α (1)

Estas expresiones evaluadas en los nudos conducen a:

uuuvvv

x yx yx y

x yx yx y

ABCABC

a

b

c

a

b

c

a a

b b

c c

a a

b b

c c

=

1 0 0 01 0 0 01 0 0 00 0 0 10 0 0 10 0 0 1

1

1

1

2

2

2

Simbólicamente: { }

( ) ( )[ ]{ }

( )q F

6 1 6 60

6 1× × ×= α (2)

De esta ecuación se obtiene: { } [ ] { }α = −F q0

1 Sustituyendo en (1)

( ){ }( )

( )[ ]( )

[ ]( )

{ }( )

( ){ }( )

( )[ ]( )

{ }( )

r x y F x y F q

r x y x y q

, ,

, ,

2 1 2 60

6 6

1

6 1

2 1 2 6 6 1

× × ×

−

×

× × ×

=

=

1 244 344

θ (3)


8

ELEMENTO FINITO RECTANGULAR

y

x

vd

cud

a

va

ua

c d

bub

v(x,y)

P(x,y)u(x,y)

vc

u c

vb

0

Nota:Los bordes delElemento Finitoson paralelos alos ejes x e y.

Siguiendo criterios similares a los empleados en el elementos finitos triangular, escogemos las siguientes funciones de comportamiento.

( )( )

u x y A B x C y D xy

v x y A B x C y D xy

,

,

= + + +

= + + +1 1 1 1

2 2 2 2

VERIFICACION DE LAS CONDICIONES DE u(x, y) y v(x, y)

− La capacidad de movimiento como cuerpo rígido se obtiene con:

A u A vB BC CD D

1 0 2 0

1 2 0

1 0 2

1 2

00

0 0

= == == − == =

;;;;

φφ

− Compatibilidad geométrica inter-elemental.

Al escoger el sistema coordenado paralelo a los bordes del elemento finito, las funciones de comportamiento pasan a ser de primer grado en dichos bordes pues en ellos x o y son constantes, lo que conduce a que los bordes se mantengan rectos en la posición deformada, satisfaciendo esta condición.

La matriz de interpolación se obtiene de forma similar a la obtenida para el Elemento Finito triangular.


9

ELEMENTO FINITO CUADRILATERAL

y

x

a

d

b

v c

u c

vd

b

vb

v a

u a

ud

b

u b

v (x,y )

P(x,y )

u(x,y )

c

0 Este elemento requiere funciones de comportamiento similares a las del elemento rectangular, pero al tener bordes inclinados con respecto a los ejes coordenados x e y no se cumple la compatibilidad inter-elemental. Para poder utilizar elementos de esta forma se emplean otros procedimientos, entre los que cabe destacar:

− Elementos finitos refinados. − Coordenadas isoparamétricas.

Materias que se estudian más adelante.

APLICACIÓN DE LOS PRINCIPIOS BASICOS DE ANALISIS ESTRUCTURAL

a) Compatibilidad Geométrica.

Las ecuaciones de la teoría de la elasticidad para pequeños desplazamientos originan las siguientes relaciones de compatibilidad geométrica.

( ) ( )

( ) ( )

( )

ε∂∂

ε ε

ε∂∂

ε ε

γ∂∂

∂∂

γ γ

x x x

y y y

xy xy

ux

x y u u x y

vy

x y v v x y

uy

vx

x y

= = =

= = =

= + =

, , , ,

, , , ,

, ,


10

Utilizando notación matricial y matrices de operadores diferenciales podemos escribir estas mismas ecuaciones en la forma:

ε

ε

γ

∂ ∂

∂ ∂∂ ∂ ∂ ∂

x

y

xy

x

yy x

uv

=

0

0

Simbólicamente: ( ){ }

( )[ ]( )

( ){ }( )

ε x y r x y, ,3 1 3 2 2 1× × ×

= ∆ (1)

Deberá notarse que la matriz [∆] contiene operadores diferenciales, por lo que no podrá estar aislada del operando sobre el cual actúa. Por otro lado tenemos: ( ){ }

( )( )[ ]( )

{ }( )

r x y x y qm m

, ,2 1 2 1× × ×

= θ (2)

m : número total de grados de libertad del elemento. (∆ : m = 6 ; � : m = 8 ; etc.)

Reemplazando (2) en (1) se tiene: ( ){ } [ ] ( )[ ]{ }ε θx y x y q, ,= ∆ (3)

b) Relaciones constitutivas

Las relaciones constitutivas para un sólido tridimensional linealmente elástico, homogéneo e isótropo y para pequeños desplazamientos, esta dada por la expresión siguiente:

( )( )

( )

γ

γ

γ

ε

ε

ε

µµ

µµ+λλλ

λµ+λλλλµ+λ

τ

τ

τ

σ

σ

σ

zx

yz

xy

zz

yy

xx

zx

yz

xy

zz

yy

xx

000000000000000000200020002


11

En que:

( )( )

( )

Lamé de Constantes

Poisson de Módulo :

=υ+

=µ

υυ−υ+

υ=λ

GE

;E

12

211

En forma matricial podemos escribir esta relación de la manera siguiente: { } [ ]{ }τ ε= c (4)

ELASTICIDAD PLANA

Al reducir el problema de elasticidad tridimensional a un problema de elasticidad plana, estamos suponiendo que una de las tres dimensiones del sólido es muy pequeña frente a las otras dos, los que nos permite efectuar simplificaciones que conducen a sólo tres componentes independientes para los vectores {τ} y {ε}. El problema plano tiene dos orígenes diferentes, conocidos como de tensiones planas o de deformaciones planas, conforme a las condiciones de tensiones o de deformaciones que existan en las dos caras paralelas separadas por el espesor del elemento.

y

z xx

e

DEFORMACIONES PLANAS

Estudiando aisladamente el elemento de espesor e, suponiendo que es representativo de toda la estructura cuyo largo es muy grande, observamos que las condiciones que le impone el resto de la estructura son las de confinar sus desplazamientos al plano (x,y), sin posibilidades de expansión según z. El problema, por consiguiente, es de deformaciones planas y las relaciones de borde en las caras adyacentes al resto de la estructura son las siguientes:


12

0

0

00

0

=τ

=τ

=γ=γ

=ε

yz

xz

yz

xz

zz

De las relaciones constitutivas se obtiene: ( )σ λ ε εzz xx yy= +

Las relaciones constitutivas para el problema de deformaciones planas se reducen a vectores de (3×1) para { τ } y { ε } y a una matriz de (3×3) para [C], pues a pesar de existir cuatro componentes no nulas de {τ}, una de ellas (σzz) no es independiente. La matriz resultante es:

( )

( )σ

σ

τ

λ µ λ

λ λ µµ

ε

ε

γ

xx

yy

xy

xx

yy

xy

=

+

+

2 0

2 00 0

con: ( )σ λ ε εzz xx yy= +

y

x

z

e

TENSIONES PLANAS

La figura, idéntica en forma a la anterior, representa un elemento que no se encuentra restringido en sus caras externas, pudiendo desplazar según z, pero sus caras estarán libres de tensiones. En consecuencia, sólo se podrán generar tensiones en el plano xy, por lo que el problema se conoce como de tensiones planas.


13

Las condiciones de borde en estas caras son las siguientes:

0

0

00

0

=γ

=γ

=τ=τ

=σ

yz

xz

yz

xz

zz

De la expresión de σzz obtenemos: ( ) ( )σ λ ε ε λ µ εzz xx yy zz= = + + +0 2

( )

ελ ε ε

λ µzzxx yy= −

+

+ 2

Lo que conduce nuevamente a tres componentes independientes para los vectores {τ} y {ε} y a una matriz constitutiva de (3×3). Particularizando la ecuación (4) para el problema plano podemos escribir la relación. ( ){ } [ ] ( ){ }τ εx y x yC, ,= (4) y substituyendo {ε(x,y)} de la ecuación (3) obtenemos: ( ){ } [ ][ ] ( )[ ]{ }τ θx y x yC q, ,= ∆ (5)


14

c) Relaciones de equilibrio: Elemento Finito triangular

x

P(x,y)

Qyc

Qxc

Qxb

QybQxa

Qy

c

b

a

y

σxx (x, y)σyy (x, y)τxy (x,y)

0

En el punto P interior del elemento finito triangular existen las tres componentes de tensión de punto. El elemento interactúa con el resto del sólido a través de fuerzas de interacción en sus puntos nodales, las que convencionalmente definimos en forma consistente con el vector {q} (esto es, equivalente a los esfuerzos en los extremos de una barra). Puesto que el elemento interactúa con el sólido solamente en sus puntos nodales, las ecuaciones de equilibrio, usando el principio de desplazamientos virtuales, las podemos escribir de la forma siguiente:

{ } { }TVE q QT

= ** : virtual

( ){ } ( ){ }TVI edxdyx yT

x yA

= ∫ ε τ*, ,


15

Con:

{ }Q

Q

Q

Q

Q

Q

Q

xa

ya

xb

yb

xc

yc

=

Igualando trabajos virtuales: TVE TVI=

{ } { } ( ){ } ( ){ }

( ){ } [ ] ( )[ ]{ }

q Q e dxdy

q

Tx y

Tx y

A

x y x y

* *, ,

* , ,*

=

=

∫ ε τ

ε θ

Y aplicando compatibilidad geométrica :

∆

(6)

Por lo tanto:

{ } { } { } [ ] ( )[ ] [ ] [ ] ( )[ ] { }q Q e q x y C x y q dxdyT T

A* * , ,=

∫ ∆ ∆θ θ

Dado que los vectores {q} y {q*} son permeables a la integral y que además {q*} es virtual y por lo tanto arbitrario, obtenemos.

{ } [ ] ( )[ ] [ ] [ ] ( )[ ][ ]

{ }Q e x y C x y dx dy

k

qT

A=

∫ ∆ ∆θ θ, ,1 244444444 344444444

(7)

[k] es la matriz de rigidez del elemento finito, la que será integrada mediante ensamble indicial a las matrices de los restantes elementos finitos de la estructura y corregida posteriormente con las condiciones de apoyo para obtener finalmente la matriz de rigidez de toda la estructura. Como se puede apreciar, una vez obtenida la matriz de rigidez de un elemento se emplean las mismas técnicas que hemos estudiado para estructuras de barra. La mayor complejidad del método reside en la búsqueda de funciones de comportamiento satisfactorias, lo cual ya ha sido logrado, disponiéndose de una amplia biblioteca de elementos finitos, de distintos tipos, para los más diversos problemas de elasticidad.


16

Los problemas algebraicos suelen ser también importantes toda vez que hay que invertir matrices con términos literales ([Fo]-1 ) y también la posterior integración. En este terreno han tenido mucho éxito los métodos numéricos de integración, entre los que destaca en forma especial la cuadratura de Gauss. Ultimamente se ha contado con un aporte notable para la solución integral de estas dificultades con el uso del programa REDUCE (Referencia: Memoria de Título Ingeniero Sr. Iván Hrepic).

VECTOR DE CARGAS

Las solicitaciones externas sobre la estructura deben ser trasladadas a los puntos nodales, para lo cual es lícito utilizar supuestos simplificatorios, tales como repartir cargas que quedan en el interior de un elemento a los nudos de éste, utilizando criterios diversos. El más utilizado de todos estos procedimientos es el que se denomina cargas consistentes y que se fundamenta en los mismos principios aplicados para obtener masas consistentes. Las cargas consistentes corresponden a cargas en los nudos del elemento que provienen de las cargas que existen en el interior del mismo elemento y que se determinan utilizando las mismas funciones de desplazamiento que las empleadas para calcular la matriz de rigidez del elemento. Para ello, supongamos un elemento finito triangular en cuyo interior actúa una carga puntual con componentes según x e y.

Qyc

Qxc

Qxb

Qyb

Qya

Qxa

c

b

a

y

x

P(x,y)Px

Py

c

b

a

y

x 0 0

Estas cargas conducen a las seis componentes nodales de la figura de la derecha, las que son estáticamente correspondientes con las primeras.


17

Usando el principio de desplazamientos virtuales y empleando las funciones de desplazamiento que sirvieron para determinar la matriz de rigidez, la equivalencia estática se puede escribir en la forma:

( ) ( ) c*cb

*ba

*ac

*cb

*ba

*a

** QyvQyvQyvQxuQxuQxuPyy,xvPxy,xu +++++=+

( )( )

( )[ ]y*

*

*

*

*

*

*

*

;,

,,

u x y

v x yx y

u

u

u

v

v

v

a

b

c

a

b

c

=

θ

En forma matricial:

( ){ } { } { }r x yPP

q QT x

y

T* *,

=

( ){ } ( )[ ]{ }r x y x y q* *, ,= θ

{ } ( )[ ]Q x yPP

T x

y=

θ ,

Distintas cargas en el mismo elemento se tratan en forma independiente, para posteriormente superponer sus efectos.


18

ELEMENTOS FINITOS REFINADOS

Para mejorar la precisión en la solución de un problema de elasticidad que se resuelve con elementos finitos, surge como primera alternativa la de refinar la malla, disminuyendo el tamaño de cada elemento, con lo que las aproximaciones muestran una convergencia a la solución exacta. Alternativamente se pueden emplear los denominados elementos finitos refinados, los que se basan exactamente en los mismos principios generales ya estudiados, pero que adicionalmente incorporan a ellos otros grados de libertad, lo que permite funciones de comportamiento de mayor grado, pudiendo obtenerse resultados satisfactorios con muchos menos elementos que los requeridos con los procedimientos convencionales. Podemos señalar que en general existen dos tipos de grados de libertad adicionales: los que provienen de la incorporación de nuevos nudos en el mismo elemento y los que provienen de la definición de grados de libertad diferentes a los traslacionales en los mismos nudos del elemento convencional. Algunos ejemplos de elementos finitos refinados son los siguientes:

Se usan polinomios completosde 2º grado.

en el centro

y

x

6 GL u6 GL v

y

x

6 GL por nudo:

4 nudos

0

0

uux

uy

vvx

vy

, , , , ,∂∂

∂∂

∂∂

∂∂

Se usan polinomios incompletos de cuarto grado, condoce términos. (Se eliminan los términos x4, x2y2, y4).

Estos ejemplos muestran los dos tipos de elementos refinados que hemos mencionado. En algunas aplicaciones específicas se suele incluir grados de libertad en nudos internos del elementos finitos, lo que aporta algunos beneficios: permite aumentar el grado del polinomio de las funciones de desplazamiento y posteriormente, una vez obtenida la matriz de rigidez del elementos finitos, estos grados de libertad internos pueden ser condensados estáticamente quedando su efecto implícito pero disminuyéndose el número de grados de libertad.


19

Un ejemplo clásico de este procedimiento lo constituye el elemento finito cuadrilateral de cinco nudos (cuatro vértices y un nudo en el centro de gravedad del cuadrilátero).

y

x 0

Para determinar la matriz de rigidez de este elemento se genera previamente las matrices de los triángulos en que queda dividido, las cuales se ensamblan formando una matriz de (10×10) que posteriormente se condensa estáticamente a (8×8) mediante la eliminación de los dos grados de libertad del nudo interior. La formulación del método para elementos refinados es idéntica a la que hemos estudiado para elementos convencionales, salvo las complicaciones algebraicas que surjan en forma específica.


20

ELEMENTOS FINITOS ESPACIALES

El problema de elasticidad espacial es estrictamente equivalente al problema plano, salvo que los desplazamientos, deformaciones y tensiones no sufren reducciones; quedando los desplazamientos con tres componentes (u,v,w) y las deformaciones y tensiones con seis componentes cada una, según se identificó en las relaciones constitutivas. A modo de ejemplo, observemos lo que ocurre con el elemento finito espacial más simple, que corresponde a la pirámide de bases triangulares.

z, w

y, v

x, u

d

c

a

b

w (x,y,z)

v (x,y,z)

u (x,y,z)

P ( x ,y,z )

( ){ }( )

( )( )( )

{ }( )

r x y z

u x y z

v x y z

w x y z

q

uuuuvvvvwwww

x

a

b

c

d

a

b

c

d

a

b

c

d

, , ;3 1 12 1×

=

=

, ,

, ,

, ,

{ }( )

{ }( )

zu

xw

;zw

yw

zv

;yv

xv

yu

;xu

;

zxzz

yzyy

xyxx

zx

yz

xy

zz

yy

xx

x

zx

yz

xy

zz

yy

xx

x

∂∂

+∂∂

=γ∂∂

=ε

∂∂

+∂∂

=γ∂∂

=ε

∂∂

+∂∂

=γ∂∂

=ε

τ

τ

τ

σ

σ

σ

=τ

γ

γ

γ

ε

ε

ε

=ε1616


21

Si: { } [ ] { }ε = ∆ r entonces:

[ ]∆ =

∂ ∂∂ ∂

∂ ∂∂ ∂ ∂ ∂

∂ ∂ ∂ ∂∂ ∂ ∂ ∂

xy

zy x

z yz x

0 00 00 0

00

0

FUNCIONES DE DESPLAZAMIENTO

( )( )( )

u x y z A B x C y D zv x y z A B x C y D zw x y z A B x C y D z

, ,, ,, ,

= + + += + + += + + +

1 1 1 1

2 2 2 2

3 3 3 3

Estas funciones permiten movimiento de cuerpo rígido por ser de primer grado y por el mismo motivo, garantizan las compatibilidades inter-elementales, pues las bases triangulares continúan siendo planas en la posición deformada. En forma general podemos señalar que las ecuaciones detalladas para los elementos finitos planos se expanden automáticamente para los problemas espaciales.


22

COORDENADAS ISOPARAMETRICAS

El desarrollo de diferentes tipos de Elementos Finitos tiene muchas limitaciones con el uso de expansiones polinómicas. Un ejemplo de esto lo constituye el elemento plano cuadrilateral, como hemos podido apreciar anteriormente. Se han desarrollado diferentes procedimientos alternativos, pero el más exitoso ha resultado ser el denominado coordenadas isoparamétricas. Consideremos, por ejemplo, un sólido tridimensional de N vértices y establezcamos la relación entre las coordenadas (x,y,z) de un punto interior cualquiera y la de esos N vértices, mediante ecuaciones del tipo siguiente:

x N x N x N xy N y N y N yz N z N z N z

n n

n n

n n

= + + += + + += + + +

1 1 2 2

1 1 2 2

1 1 2 2

LLL

En que: ( )N N= ξ ψ, ,η tal que:

( )Ni ji ji j j jξ η ψ, ,

;;

=≠=

01

Como se puede apreciar, (ξ , η , ψ) representan un nuevo sistema coordenado que pasa a substituir a (x, y, z), que con la condición indicada originan identidades al evaluarse las coordenadas en los vértices del sólido.


23

Por ejemplo, supongamos un problema plano para describir un cuadrilátero.

ξ=1 ξ= -1

η

ξ

η = 1

η = -1

y

x

3

1

2

4

β/2 β/2

α /2

α /2

0

Los ejes ξ y η son variables entre los valores -1 y 1 y se construyen mediante la rotación relativa de un eje que pivotea en torno a la intersección de dos lados opuestos. Con este sistema de variables los bordes pueden tener la expansión: x A A A Ay A A A A

= + + += + + +

1 2 3 4

5 6 7 8

ξ η ξ ηξ η ξ η

(1)

Que para ξ o η constante (caso de los bordes) representa ecuaciones de primer grado consistentes con bordes rectos. Si agregamos la condición que esta función satisfaga las coordenadas de los cuatro vértices, habremos logrado una descripción geométrica exacta de la figura. Para la expansión x, por ejemplo, tendríamos lo siguiente: x A A A Ax A A A Ax A A A Ax A A A A

1 1 2 3 4

2 1 2 3 4

3 1 2 3 4

4 1 2 3 4

= − − += + − −= + + += − + −


24

Determinando A1, A2, A3 y A4 en términos de las coordenadas x1, x2, x3 y x4, sustituyendo en la ecuación (1) y factorizando, se obtiene la siguiente expresión.

( )( ) ( )( ) ( )( ) ( )( ) 441

341

241

141

4321

11111111 xxxxxNNNN 44 844 7644 844 7644 844 7644 844 76

η+ξ−+η+ξ++η−ξ++η−ξ−= Notar que en los vértices se cumple:

Vértices ξ η x yx yx yx yx y

1 1 12 1 13 1 14 1 1

1 1

2 2

3 3

4 4

− −+ −+ +− +

Si evaluamos Ni en los vértices j, observaremos que se cumple exactamente la condición impuesta anteriormente, es decir:

( )Ni ji ji j jξ η,

;;

=≠=

01

El método es aun más interesante si se piensa que no hay restricción alguna para considerar bordes curvos. Por ejemplo, veamos el caso de un cuadrilátero genérico de lados parabólicos:

η = 1η

7 4

3

ξ = 1

ξ = -1 ξ8

6

5

η = -1

1 2

En este cuadrilátero generalizado se requieren tres puntos en cada borde para describir una parábola. Los ejes ξ y η serán curvos, acomodándose con sus valores -1 y 1 a los bordes que corresponden a su generación.


25

La expanción para las coordenadas requiere ocho términos y que sea una ecuación de segundo grado en ξ o η para los bordes. El polinomio de segundo grado tiene seis términos lo que nos conduce a incorporar dos términos de tercer grado que deriven en segundo grado para ξ y η constantes, lo que origina que: x A A A A A A A A= + + + + + + +1 2 3 4 5

26

27

28

2ξ η ξη ξ η ξ η ξη El resto del procedimiento es similar al caso de bordes rectos. Si se deseara desarrollar un cuadrilátero de bordes curvos que correspondan a parábolas cúbicas, se sigue exactamente el mismo procedimiento, eligiendo nudos en los tercios de los bordes, totalizando doce nudos para el elemento y se emplean expansiones completas de tercer grado (diez términos) más dos términos de cuarto grado que derivan en tercer grado para ξ o η constantes y que por lo tanto son iguales a ξ3η y ξη3. El resto del procedimiento es análogo.

FORMULACION DE LAS ECUACIONES CINEMATICAS.

Con el mismo criterio con que hemos descrito la geometría del cuerpo, podríamos describir sus desplazamientos en función de los desplazamientos nodales, es decir:

u M u M u M uv M v M v M vw M w M w M w

n n

n n

n n

= + + += + + += + + +

1 1 2 2

1 1 2 2

1 1 2 2

KKK

en que: ( )M M= ξ η ψ, , con:

con: ( )

=≠

=ψηξji;ji;

,,M jjji 10


26

Esta ecuación escrita en forma matricial conduce a:

( )[ ]

ψηξθ

=

n

n

n

n

n

n

w

wwv

v

vu

u

u

,,

MMMMMM

MMM

wvu

M

M

M

4444444444444 34444444444444 21LLLLLLLLL

2

1

2

1

2

1

21

21

21

000000000000000000

Al identificar los términos observamos que esta ecuación corresponde a: ( ){ } ( )[ ]{ }r x y z q, . , ,= θ ξ η ψ

Lo que indica que al determinar las funciones Mi , automáticamente obtenemos la matriz de interpolación. Los desplazamientos cumplirán las condiciones de movimiento de cuerpo rígido si en sus expansiones poseen al menos el desarrollo completo de los términos de primer grado. Consideremos, por ejemplo, los desplazamientos u y determinemos los requerimientos que originan las expansiones mínimas de primer grado. u A Bx Cy Dz= + + + Por otro lado: u M u M u M un n= 1 1 2 2+ + +K Por lo tanto:

( ) ( )( ) ( )

u M A Bx Cy Dz M A Bx Cy Dz

M A Bx Cy Dz M A Bx Cy Dz

= 1 1 1 1 2 2 2 2

3 3 3 3 4 4 4 4

+ + + + + + +

+ + + + + + + +


27

Factorizando: u A M B M x C M y D M zi i i i i i i= ∑ ∑ ∑ ∑+ + +

Ecuación que comparada con el desarrollo en primer grado resulta en:

M

M x x

M y y

M z z

i

i i

i i

i i

∑∑∑∑

=

=

=

=

1

Las tres últimas ecuaciones conducen a identificar exactamente las funciones Mi con las Ni . Estas funciones Ni que describen la geometría podrán ser utilizadas para describir los desplazamientos si su sumatoria es igual a uno. La compatibilidad geométrica inter-elemental quedará automáticamente satisfecha pues en la evaluación de las funciones Ni hemos forzado el cumplimiento de la geometría de los bordes. El hecho de poder emplear las mismas funciones para describir geometría y cinemática de un sólido conduce al nombre de coordenadas isoparamétricas. Observemos el ejemplo del cuadrilátero de bordes rectos que hemos desarrollado con anterioridad y veamos a qué conduce la sumatoria de los Ni.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

( )

N

N

N

i

i

i

∑

∑∑

= − − + + − + + + + − +

= + + +

=

14

1 114

1 114

1 114

1 1

14

1 1 1 1

1

ξ η ξ η ξ η ξ η

Lo que demuestra que estas funciones son isoparamétricas. En la determinación de la matriz de rigidez del elemento finito se tiene una expresión del tipo:

[ ] [ ] [ ][ ] [ ] [ ] [ ][ ]K C dVT

v= ∫ ∆ ∆θ θ

Ecuación que debe ser referida a un mismo sistema coordenado para poder desarrollar la integración. Los límites de la integración surgen de inmediato de la descripción geométrica del sólido; la matriz de interpolación ya se encuentra transformada a las nuevas


28

coordenadas y la matriz constitutiva no depende del sistema coordenado. Por lo tanto sólo requerimos transformar la matriz de operadores diferenciales y el diferencial de volumen. Del cálculo diferencial se obtienen las siguientes expresiones:

{ } [ ]{ } [ ]

{ } [ ] { }

∂ ∂ξ∂ ∂η∂ ∂ψ

∂ ∂ξ ∂ ∂ξ ∂ ∂ξ∂ ∂η ∂ ∂η ∂ ∂η∂ ∂ψ ∂ ∂ψ ∂ ∂ψ

∂ ∂∂ ∂∂ ∂

∂ ∂

∂ ∂

=

→ =

=−

x y zx y zx y z

xyz

J J

J

*

*

; : Jacobiano

1

{ } [ ] { }*Jzy

x

x ∂><=∂><=

∂∂∂∂

∂∂

><=∂∂ −1001001001

{ } [ ] { }*Jzy

x

y ∂><=∂><=

∂∂∂∂

∂∂

><=∂∂ −1010010010

{ } [ ] { }*Jzy

x

z ∂><=∂><=

∂∂∂∂

∂∂

><=∂∂ −1100100100

Ecuaciones que nos permiten expresar cada una de los componentes de la matriz [∆] en función del nuevo sistema coordenado. Por su parte, el diferencial de volumen es igual a:

[ ]dV dx dy dz Det J d d d= = ξ η ψ


29

COORDENADAS NATURALES

ELEMENTOS TRIANGULARES

P(x,y)

L1=0 L2=0

L3=0

3 y

2

1

x 0

Como se recordará, las coordenadas isoparamétricas son aquellas que junto con describir la geometría del cuerpo permiten también describir sus desplazamientos, lo que requiere que la suma de Ni para i de 1 a n, sea igual a uno.

Nii

n

=∑ =

1

1

Esta condición debe ser verificada una vez determinados los Ni , salvo que se escojan funciones que a priori la satisfagan, como es el caso de las denominadas coordenadas naturales (también conocidas como coordenadas de área) para elementos de forma triangular. En la figura anteior observamos un triángulo en cuyo interior identificamos un punto P(x,y). Este punto puede quedar también descrito mediante nuevas coordenadas L1, L2 y L3, que se definen de la manera siguiente:

La Pa

La Pa

La Pa

a1 2 323

12331

12312

123= = = =

``

``

``

`∆∆

∆∆

∆∆

; ; ; area

De donde se desprende que: L L L1 2 3 1+ + = Lo que permite señalar que estas nuevas coordenadas son isoparamétricas.


30

Los bordes del elemento triangular resultan con ecuaciones referidas a este nuevo sistema coordenado, que se pueden escribir de la siguiente forma:

Borde EcuaciónLLL

3

1

2

1 2 02 3 03 1 0

− =− =− =

Los nudos tienen las siguientes coordenadas:

Nudo x y L L Lx yx yx y

1 2 3

1 1

2 2

3 3

1 1 0 02 0 1 03 0 0 1

Si escribimos la ecuación que define la geometría del elemento en la forma: x A L A L A L= + +1 1 2 2 3 3 y satisfacemos el paso por los tres vértices, obtenemos:

x Ax Ax A

1 1

2 2

3 3

===

por lo tanto: x L x L x L x= + +1 1 2 2 3 3 N L N L N L1 1 2 2 3 3= = =; ;

La matriz de interpolación será la que establece la siguiente relación:

uv

L L LL L L

uuuvvv

=

1 2 3

1 2 3

1

2

3

1

2

3

0 0 00 0 0


31

Esta formulación permite desarrollar elementos refinados triangulares en forma relativamente simple, si se le compara con iguales procedimientos en el sistema cartesiano original.

EJEMPLOS DE APLICACION

ELEMENTO TRIANGULAR DE SEIS NUDOS

L 1=0 L 2=0

L 3=0

3 y

5

2

4

6

1

0 x

Los nudos 4, 5 y 6 se ubican en los puntos medios de los lados. La definición de un tercer punto en cada lado del triángulo exige la definición de funciones de segundo grado cuyo coeficientes se ajustan de tal manera que al actuar sobre las coordenadas de los nudos deriven en ecuaciones de primer grado, pero que al actuar sobre los desplazamientos nodales permitan que los lados describan una curva de segundo grado en la posición deformada. La geometría queda descrita en la forma: x A L A L A L A L L A L L A L L= + + + + +1 1 2 2 3 3 4 1 2 5 2 3 6 3 1

Nudo x y L L Lx yx yx yx yx yx y

1 2 3

1 1

2 2

3 3

4 41

21

2

5 51

21

2

6 61

21

2

1 1 0 02 0 1 03 0 0 14 05 06 0


32

Obligando a que esta función pase por los seis nudos, se tiene:

x Ax Ax A

x A A A

x A A A

x A A A

1 1

2 2

3 3

4 1 2 4

5 2 3 5

6 1 3 6

12

12

14

12

12

14

12

12

14

===

= + +

= + +

= + +

Resolviendo el sistema de ecuaciones, reemplazando y expresando: L L L3 1 21= − − se tiene: x N x N x N x N x N x N x= + + + + +1 1 2 2 3 3 4 4 5 5 6 6 en que:

( )( )

( )( )

( )( )

N L L

N L L

N L L L L

N L L

N L L L

N L L L

1 1 1

2 2 2

3 1 2 1 2

4 1 2

5 2 1 2

6 1 1 2

2 1

2 1

1 1 2 2

4

4 1

4 1

= −

= −

= − − − −

=

= − −

= − −

Pudiendo demostrarse que:

N ii=∑ =

1

6

1


33

ELEMENTO TRIANGULAR DE CUATRO NUDOS.

x

L1=0 L2=0

L3=0

3 y

4

2

1

0 Este elemento se resuelve con las mismas técnicas anteriores, pero con la siguiente ecuación de descripción geométrica:

44332211315324332211 xNxNxNxNxLLALLALALALAx +++=++++= La que origina ecuaciones de primer grado para los bordes 12 y 13 y de segundo grado para el borde 23. Resolviendo en forma similar al caso anterior se obtiene:

( )( )

( )214

33

122

211

4

21

21

LLN

LN

LLN

LLN

=

=

=

=

−

−

ELEMENTO TRIANGULAR DE CINCO NUDOS.

3

1

0

L1=0 L2=0

L3=0

y

x

5 2

4

x A L A L A L A L L A L L= + + + +1 1 2 2 3 3 4 2 3 5 1 3

5544332211 xNxNxNxNxNx ++++= ( )( )( )

315

324

2133

322

311

4

4

221

21

21

LLN

LLN

LLLN

LLN

LLN

=

=

=

=

=

−−

−

−


34

PIRAMIDE DE BASE TRIANGULAR

Las coordenadas naturales o de área se expanden con facilidad al caso de una pirámide de base triangular, la que requiere cuatro términos (L1, L2, L3 y L4) que miden el cuociente entre los volúmenes de las cuatro pirámides parciales en que queda dividida la pirámide original al definir un punto en su interior y unirlo con los cuatro vértices y el volumen de la pirámide total.

LVPV

LVPV

LVPV

LVPV

1

2

3

4

2341234

1341234

1241234

1231234

=

=

=

=

L L L L1 2 3 4 1+ + + =

P

4

3

2

1


35

COMENTARIOS

Los temas estudiados cumplen el propósito de introducción al Método de Elementos Finitos, brindando alguna capacidad de solución de problemas. Para entrar en más detalle se podrá consultar la bibliografía especializada que se adjunta. En todo caso, conviene resaltar algunos aspectos prácticos importantes. El primero de ellos se refiere a la determinación de tensiones internas, las que deberán ser calculadas en forma consistente con las simplificaciones del método, no pudiendo aceptarse como válidas aquellas que aparentemente se podrían deducir de la ecuación.

( ){ } [ ] ( ){ } [ ][ ] ( )[ ]{ }τ ε θx y x y x yC C q, , ,= = ∆ En virtud que el equilibrio sólo está garantizado en los nudos del elemento finito. Por este motivo, la práctica habitual sugiere que una vez conocidos los desplazamientos nodales se evalúen las tensiones en esos mismos nudos, promediándolas posteriormente y estableciendo que dicho promedio es válido para todo el elemento como un valor puntual en su centro de gravedad. Para obtener el estado de tensiones en una sección cualquiera del cuerpo, se deberá pasar una curva, con buen criterio, por tales centros de gravedad. El segundo aspecto de importancia se refiere al hecho que algunas singularidades no quedan bien resueltas por el método. Por ejemplo, debido al mismo argumento anterior, los bordes libres no provocarían tensiones nulas y éstas disminuirían su valor en la medida que los elementos sean más pequeños. Esto exige muchas veces utilizar refinamientos locales de las mallas de elementos finitos, para poder lograr soluciones prácticas satisfactorias.


36

BIBLIOGRAFIA

1.- The Finite Element Method in Engineering Science. C. ZIENKIEWICZ 2.- Theory of Matrix Structural Analysis. S. PRZEMIENIECKI 3.- Finite Element Method. Editor IVAR HOLAND 4.- Numerical Methods in Finite Element Analysis. KLAUS-JÜRGEN BATHE 5.- The Finite Element Technique. A. BREBBIA 6.- El Método de Elementos Finitos en Ingeniería Civil. Publicación Sección Estructuras Dpto. Ingeniería Civil U. de Chile. T. GUENDELMAN, M. KATANELLA, R. SARAGONI, M. SARRAZIN, P. VIVERO

Documents

Capítulo 12: El Método de Elementos Finitospolisoft.cl/clubpenrose/data/tg/tg_metodoelementosfinitos.pdf · El método de elementos finitos es una técnica de discretización física