Alberto Cardona, Víctor Fachinotti Cimec-Intec (UNL/Conicet), Santa Fe, Argentina Introducción al...

Alberto Cardona, Víctor FachinottiCimec-Intec (UNL/Conicet), Santa Fe, Argentina

Introducción al Método de los Elementos Finitos


2

Introducción

Modelos matemáticos en ciencia e ingeniería Ecuaciones algebraicas, diferenciales o integrales

El desarrollo de las computadoras permitió usar estos modelos para resolver problemas prácticos. Se pueden simular y resolver sistemas altamente complicados en ciencia e ingeniería.

Permiten:

1. Reducir la necesidad de experiencias con modelos y prototipos (caras y lentas).

2. Comparar fácilmente distintas alternativas de diseño para llegar al óptimo ingenieril.

Disciplinas relacionadas:

1. CAD: Computer Aided Design

2. CAE: Computer Aided Engineering

3. CAM: Computer Aided Manufacturing


3

Introducción (cont)

Método de los Elementos Finitos (MEF) : técnica general para hallar soluciones numéricas de sistemas de ecuaciones diferenciales e integrales.

Origen: ingeniería estructural, años 50/60, para solución de ecuaciones diferenciales en derivadas parciales en elasticidad.

Su aplicación se generalizó, integrado a sistemas de CAD/CAE.

Modelos matemáticos

Soluciones analíticas

Soluciones numéricas

Sólo casos simples

Casos prácticos, en computadora


4

Introducción (cont)

Aplicaciones del MEF:

• Ingeniería estructural

• Resistencia de materiales

• Mecánica de fluidos

• Ingeniería nuclear

• Electromagnetismo

• Campos eléctricos

• Propagación de ondas

• Conducción del calor

• Procesos de convección – difusión

• Ingeniería de petróleo

• Procesos de reacción – difusión

• . . . . . . . .


5

Método de diferencias finitas vs. MEF

•Idea básica de un método numérico para resolver ecuaciones diferenciales:

•El problema discreto puede resolverse en una computadora !!!

•Diferencias finitas: método numérico clásico para resolver ecuaciones diferenciales donde el problema discreto se obtiene reemplazando:

Continuo con infinitos GDL

Problema discreto c/ # GDL finito

Discretización

Solución discreta, aproximada, con N GDL

1 nnu uux x


6

Discretización en el MEF

1. Reformulación de la ecuación diferencial en un problema variacional equivalente.

Ej: en ecuaciones elípticas, en casos simples, toma la forma de problema de minimización:

donde:

V: conjunto de funciones admisibles

es un funcional

las funciones a menudo representan cantidades que varían en forma continua (ej.. desplazamiento de un cuerpo elástico, temperatura, etc.)

F(v) es la energía total asociada a v.

(M) es una caracterización equivalente de la solución de la ecuación diferencial como aquella función en V que minimiza la energía total del sistema considerado.

(M) Hallar / ( ) ( )u V F u F v v V

:F V ( )F v v V v V


7

Discretización en el MEF (cont)

2. En general, la dimensión de V es infinita (las funciones de V no pueden expresarse a través de un número finito de parámetros). Luego, (M) no puede resolverse en forma analítica.

Para hallar una solución, la idea del MEF es reemplazar por un conjunto de funciones simples que dependen de un número finito de parámetros:

Este problema es equivalente a un sistema de ecs. algebraicas. Se espera que sea una aproximación suficientemente buena de , la solución de (M).

3. Usualmente elegimos:

En este caso, corresponde al método clásico de Ritz-Galerkin

V hV

(M ) Hallar / ( ) ( )h h h h hu V F u F v v V

huu

/h h h hV V v V v V

(M )h


8

Discretización en el MEF (cont)

4. Característica particular del MEF: funciones de son funciones polinomiales por tramos

hV

Veremos más adelante que se pueden hacer formulaciones variacionales más

generales que el problema de minimización; ej.: métodos de Galerkin

Pasos para resolver por MEF:

1. Formulación variacional del problema

2. Discretización por MEF: construcción del espacio dimensional finito

3. Solución del problema discreto

4. Implementación del método en computadora: programación

hV


9

Comentarios

• Existen distintas formulaciones variacionales dependiendo, por ejemplo, de

la elección de las variables independientes

• La elección del subespacio de aproximación de dimensión finita , o sea la

elección del tipo de elemento finito, está influenciada por:

– Formulación variacional

– Requerimientos de precisión

– Propiedades de regularidad de la solución exacta

– . . .

• Para resolver el problema discreto necesitamos algoritmos de optimización

y/o algoritmos de solución de grandes sistemas de ecs. algebraicas lineales o

no lineales

hV


10

Ventajas del MEF respecto de DF

– Tratamiento de geometrías complicadas

– Condiciones de borde generales

– Propiedades materiales no lineales o variables

– La estructura clara del método permite crear códigos multipropósito

generales

– Fundamentos teóricos sólidos

– Confiabilidad

– Posibilidad de estimación de error


11

1. Introducción al MEF para problemas elípticos

• Problemas elípticos “modelo” y solución por MEF

– Problema simple 1-D

– Generalización a 2-D

• Propiedades básicas del método


12

(1.1) Formulación variacional del problema 1-D

• Sea el problema de valores de frontera:

donde es una función continua dada.

• Integrando 2 veces, vemos que el problema tiene solución única u.

• (D) puede describir cualquiera de los siguientes problemas de Mecánica del continuo:

A. Barra elástica

B. Cuerda elástica

C. Conducción del calor en una barra

'' ( ) , 0 1(D)

(0) (1) 0

u f x x

u u

' ; ( )du

u f xdx


13

A) Barra elástica

Barra elástica sujeta en ambos extremos y sometida a carga axial de intensidad f(x)

Bajo hipótesis de pequeños desplazamientos y material elástico lineal:

donde E es el módulo de elasticidad. Asumiendo E=1

' Ley de Hooke

' Ec. de equilibrio

(0) (1) 0 Cond. de borde

Eu

f

u u

'' ( ) , 0 1(D)

(0) (1) 0

u f x x

u u


14

B) Cuerda elástica

Cuerda elástica sujeta en ambos extremos, con tensión unitaria y sometida a carga transversal de intensidad f(x)

Bajo hipótesis de pequeños desplazamientos y material elástico lineal:

'' ( ) , 0 1(D)

(0) (1) 0

u f x x

u u


15

C) Conducción de calor en una barra

Barra sometida a fuente de calor distribuida f(x), con temperatura nula en ambos extremos.

Bajo condiciones estacionarias y material lineal:

donde k es el conductividad. Asumiendo k=1

' Ley de Fourier

' Ec de equilibrio

(0) (1) 0 Cond de borde

q ku

q f

u u

'' ( ) , 0 1(D)

(0) (1) 0

u f x x

u u


16

Problemas de minimización y variacional

Veremos que la solución al problema (D) es también solución del problema de minimización (M) y de un problema variacional (V)

Para formular (M) y (V) introducimos nueva notación.

1. Producto interno (v,w):

para funciones reales acotadas continuas por tramos v, w.

2. Espacio lineal V:

3. Funcional lineal

1

0( , ) ( ) ( )v w v x w x dx

/ funcion continua sobre 0,1 ,

' es continua p/tramos y acotada en 0,1 ,

(0) (1) 0

V v v

v

v v

:F V 1

( ) ( ', ') ( , )2

F v v v f v


17

Problemas de minimización y variacional (cont)

Problema (M):

Problema (V):

Nota: en el contexto de los problemas (A) y (B)es la energía potencial total asociada al desplazamiento v

es la energía elástica interna

es el potencial de cargas

Problema (M): principio de mínima energía potencial en Mecánica

Problema (V): principio de trabajos virtuales en Mecánica

Veremos la equivalencia de los problemas (D), (V) y (M)


(V) Hallar / ( ', ') ( , )u V u v f v v V

1( ', ')

2v v

( , )f v

( )F v


18

La solución de (D) es solución de (V)

1. Multiplicamos por una función arbitraria (func. de test).

Integramos sobre (0,1):

2. Integramos por partes el lado izquierdo, y usamos v(0)=v(1)=0:

3. Al ser v arbitraria:

o sea, u es solución de (V)

''u f v V

( '', ) ( , )u v f v

( '', ) '(1) (1) '(0) (0) ( ', ') ( ', ')u v u v u v u v u v

( ', ') ( , )u v f v v V (1.1)


19

Los problemas (M) y (V) tienen la misma solución

1. Sea u solución de (V). Sea

Luego:

o sea, u es solución del problema (M).

• Sea u solución de (M). Luego

Definiendo:

Por (*) tiene un mínimo en . Luego :

o sea, u es solución del problema (V).

y / .v V w v u v u w w V

( ) 0 por (1.1)

0

1( ) ( ) ' ', ' ' ( , )

21 1

', ' ( , ) ( ', ') ( , ) ', ' ( )2 2

F u

F v F u w u w u w f u w

u u f u u w f w w w F u

( y : ) ( (*) )F u F u vv V

21( ) ( ) ( ', ') ( ', ') ( ', ') ( , ) ( , )

2 2g F u v u u u v v v f u f v

( )g 0

( ) min en 0

'(0) ( ', ') ( , ) 0g

g u v f v


20

La solución de (V) es única

1. Sean

2. Sustrayendo, y eligiendo

lo cual muestra que:

3. Usando la condición de borde:

o sea, la solución a (V) es única.

11 2 2, soluciones de (V ,): y u u Vu u

1

2

( ' , ') ( , )

( ' , ') ( , )

u v f v v V

u v f v v V

1 2v u u V 1 21 20

( ' ' ) 0u u dx

1 2 1 2' ( ) ' ( ) ( ) ' ( ) 0 (0,1)u x u x u u x x

1 2 1 2(0) (0) 0 ( ) ( ) 0,1u u u x u x x


21

Equivalencia de soluciones a (D), (V) y (M)

1. Hasta ahora hemos visto que si u es solución a (D), luego es solución a los problemas equivalentes (V) y (M) :

Mostraremos que si u es solución de (V), luego u satisface (D).

2. Sea

Asumimos existe y es continua. Integ.p/partes y usando

Como es continua, luego:

o sea, u es solución de (D).

( ) ( ) ( )D V M

1 1

0 0/ ' ' 0u V u v dx fv dx v V

''u (0) (1) 0v v 1 1 1

00 0

1

0

'' ' 0

( '' ) 0

u v dx fv dx u v

u f v dx v V

( '' )u f

( '' ) ( ) 0 0 1u f x x


22

Equivalencia de soluciones (D), (V) y (M)

Resumiendo: Hemos demostrado que

1. La solución de la ecuación diferencial es solución de un problema

variacional

2. La solución del problema variacional es también solución de un problema de

minimización y viceversa

3. La solución del problema variacional es única

4. Si se cumple un requisito de regularidad (u” continua), la solución del

problema variacional es también solución de la ecuación diferencial

Notar: Las soluciones a los problemas variacional y de minimización vistos hasta ahora tienen dimensión infinita, no pueden hallarse en computadora.

Veremos ahora cómo el MEF construye aproximaciones de dimensión finita a las soluciones de (V) y (M).


23

(1.2) MEF p/problema modelo c/funciones lineales p/tramos

1. Construiremos un subespacio de dimensión finita del espacio , consistente en funciones lineales p/tramos.

2. Sea

una partición del intervalo (0,1) en subintervalos de longitud

La cantidad es una medida de la densidad de la partición.

hV V

0 1 10 1M Mx x x x

1( , )j j jI x x

1, 1, 1j j jh x x j M max jh h


24

Subespacio de funciones lineales por tramos

1. Sea

2. Para describir elegimos los valores

/ es lineal en cada subintervalo

es continua en el 0,1

(0) (1) 0

h jV v v I

v

v v

Ejemplo de hv V

Notar que hV V

( ) en los nodos , 1,j j jv x x j M

hv V


25


3. Definimos funciones de base

4. Toda función puede ser escrita en forma única como combinación lineal de las funciones de base :

Luego, es un espacio vectorial lineal de dimensión M con base:

( ) , 1,j hx V j M

1 si ( ) , 1,

0 si j i ij

i jx i j M

i j

: función continua lineal por tramos que verifica

la propiedad delta.

( )j x

hv V( )i x

1

( ) ( ) 0,1 , ( )M

i i i ii

v x x x v x

hV 1

M

i i


26

MEF para problema modelo (D)

1. Formulación como problema de minimización discreto:

2. De la manera ya vista, es equivalente al problema variacional discreto:

3. Notar, que si satisface (1.2), luego en particular

(además, si se cumple (1.3), luego vale (1.2) )

4. Siendo , escribimos (1.3) en la forma:

Sistema de M ecuaciones algebraicas lineales con M incógnitas

h(M ) Hallar / ( ) ( )h h h hu V F u F v v V

h(V ) Hallar / ( ', ') ( , )h h h hu V u v f v v V

(método Ritz)

(método Galerkin)(1.2)

( ', ') ( , ) 1,h j ju f j M h hu V

(1.3)

hv V

1

( ) ( ) , ( )M

h i i i h ii

u x x u x

1

', ' , 1, ,M

i i j ji

f j M

i


27

Los elementos pueden calcularse fácilmente.

Notar:

Luego A es tridiagonal

Forma matricial

A b1 1 1 2 1 1 1

2 1 2 2 2 2

1

( ', ') ( ', ') ( ', ')

( ', ') ( ', ')

( ', ') ( ', ')

M

M M M M M

b

b

b

A : matriz de rigidezb : vector de cargas

', 'ij i ja

( , )i ib f

', ' 0 si 1i j i k

(1.5)


28

Sistema de ecuaciones

Entonces:

1. A es simétrica pues:

2. A es definida positiva pues siendo luego:

Como

Entonces:

1

1

1

2 21 1

( 1) 1 2

1 1 1 1( ', ') 1,

1 1( ', ') 2,

j j

j j

j

j

x x

jj j j x xj j j j

x

j j j j xj j

a dx dx j Mh h h h

a dx j Mh h

', ' ', 'i j j i

1

( ) ( )M

i ii

v x x

, 1 1 1

0 casi siempre ', ' ', ' ', ' 0

0 solo si ' 0

M M M

i i j j i i j ji j i j

v vv

(0) 0 luego ( ', ') 0 0 o sea: 0 2,jv v v v j M

0 , 0M Tη Aη η ηA simétrica y definida positiva


29

Propiedades sistema de ecuaciones (1.5)

1. A sim > 0 A es no singular y el sistema (1.5) tiene solución única

2. A es rala, o sea pocos elementos de A son distintos de cero.

Esto se debe a que tiene soporte local ( en un intervalo pequeño,

interfiriendo con pocas funciones ).

3. Si la partición es uniforme, y logramos

j 0j k

1

1jh hM

2 1

1 2 11

1

2 1

1 2

h

0

A

0


30

(1.3) Estimación de error del MEF para problema modelo

Sean

Recordando que

y como

donde es el error de la aproximación

Norma asociada al producto escalar ( , ):

Desigualdad de Cauchy:

'' ( ) , 0 1

(0) (1) 0 solucion de (D)

u f xu

x

u u

h Hallar solucion de (V / ( ', ') ( )) ,h h h hh u V u v f v vu V

( ', ') ( , )u v f v v V

hV V ', ' 0hu u v v V

hu u

Ecuación del error

1

11 222

0,w w w w dx

,v w v w

Veremos que en cierto modo, es la mejor aprox posible a la sol exacta uhu


31

Estimación de error del MEF para problema modelo

Teo 1.1)

D) Sean . Reemplazando v por w en la ecuación del error:

Luego

Para lograr una estimación cuantitativa del error, usamos una elegida convenientemente. Por el resultado anterior:

' 'h hu u u v v V

arbitraria y h h hv V w u v V

ecuacion del error

Cauch

2

y

0

' ', ' ', ' ', '

', ' ' '

h h h h h h

h h

u u u u u u u u w u u u u w

u u u v u u u v v V

' 'h hu u u v v V

h hu V

' 'h hu u u u


32

Estimación del error

Haremos interpolante de , o sea:

En Análisis Numérico, se ve que:

Luego, por el teorema y (1.12):

Mediante análisis detallado, se puede mostrar :

Notar:• No necesitamos construir explícitamente , sino sólo la estimación de error del

interpolante.• es una deformación o tensión, y tiene interés práctico la estimación de su error

h hu V u

( ) ( ) 0, 1h j ju x u x j M

2

0 10 1

( ) ( ) max ( ) ( ) ( ) max ( )8h h y

y

hu x u x h u y u x u x u y

0 1max ( )h

yu u h u y

2hu u o h

hu

u


33

(1.4) MEF para la ecuación de Poisson

Sea el problema:

con

Puede corresponder a muchos problemas físicos: conducción de calor, potencial electromagnético, desplazamiento de una membrana fija en la frontera, etc.

Teorema de la divergencia:

con:

Fórmula de Green:

en (D)

0 sobre

u f

u

2 2

2 21 2

u uu

x x

div d ds

A x A n

1 11 2

2 21 2

divA nA A

A nx x

A A n

v w d v w ds v w d x n x


34

Formulación variacional

1. Si u satisface (D), luego es solución de:

con

D) Multiplicamos (D) por una arbitraria e integramos

2. Como en el caso 1-D, si u suficientemente regular :

con

Energía potencial total:

(V) Hallar / ( , ) ( , )u V a u v f v v V

( , ) (forma bilineal) ( , )a u v u v d f v f v d

x x

1 2

: (i) es continua en ; (ii) , continuas a trozos en ; v v

V v vx x

hv V(iii) 0 sobre v

0 ( 0 en ) ( , )

( , ) ( , )

v a u v

f v v u d v u ds u v d a u v

x n x

( ) ( ) ( )D V M


1( ) ( , ) ( , )

2F v a v v f v


35


1. Construiremos un subespacio de dimensión finita consistente en funciones lineales p/tramos.

2. Asumimos poligonal. Sea una triangulación de : con

triángulos no superpuestos

y de forma que no haya ningún vértice de algún triángulo ubicado sobre el

lado de otro.

El parámetro de malla mide la densidad de triangulación.

hV V

1 2, ,h mT K K K /iK

1 2

h

mK T

K K K K

hK T

K

max diam( )hK T

h K


36

Subespacio de funciones lineales por tramos (cont)

3. Definimos

Notar

4. Parámetros p/def.

Excluimos los nodos de frontera pues

5. Def. funciones de base :

: (i) es continua en ; (ii) es lineal para ; h hKV v v v K T

(iii) 0 sobre v

hV V

: valores ( ) de en los nodos 1, de h i i hv V v N v N i M T

0 sobre .v

1 si ( ) , 1,

0 si j i ij

i jN i j M

i j

( ) /j hV x

Soporte de ( ) / con nodo j h jT N x x x


37

Subespacio de funciones lineales por tramos (cont)

6. Toda función tiene luego la representación:

7. Formulamos luego el siguiente MEF p/el problema de Poisson (D):

8. De manera similar al caso 1-D, este problema es equivalente a resolver el sistema de ecuaciones:

donde ahora:

9. Nuevamente A es simétrica, definida positiva y no singular, con lo cual el sistema de ecuaciones tiene solución única. Además, A es rala pues:

hv V

1

( ) ( ) , ( )M

i i i ii

v v N

x x x

h(V ) Hallar / ( , ) ( , )h h h hu V a u v f v v V

A b

, ( , )ij i j i j i i ia a d b f f d

x x

( )i h iu N

0 si , al mismo triánguloij i ja N N


38

Noción de mejor aproximación

1. es la mejor aproximación a la solución u en el sentido:

con

2. En particular, si encontramos tal que

podremos probar convergencia. Ejemplo, usando el interpolante

tendremos: con C>0 cte independiente de h, que depende de:• tamaño derivadas segundas de u• menor ángulo de los triángulos

3. Puede mostrarse luego:

4. Así, si u es sufic. regular, el error y su gradiente tienden a cero con

h hu V

h hu u u v v V

1

1 222,v a v v v d

x

h hu V

h hu u u u

( ) ( ) 1,h j ju N u N j M

hu u Ch

hK T

1

2 2 2h hu u u u d Ch

x

0h


39

Cálculo matriz rigidez

1. Los elementos son calculadas por suma de contribuciones de los distintos triángulos:

Notar que

2. Sean los nodos del triángulo K. Luego, la matriz de rigidez del elemento K:

3. La matriz de rigidez global A es armada luego en 2 etapas:1. Cálculo de las matrices de rigidez elementales

2. Sumatoria de las contribuciones de cada elemento (ensamble)

El vector de cargas b es armado de la misma manera.

,ij i ja a

, ,h h

i j i j i j K i jKK T K T

a d d a

x x

, 0 sólo si ,K i j i ja N N K

, y i j kN N N

, , ,

, ,

sim. ,

K i i K i j K i k

K K j j K j k

K k k

a a a

a a

a

A


40

Cálculo matriz rigidez elemental

1. Trabajamos con las restricciones de las funciones de base al triángulo K:

2. es una función lineal /

3. Si es una función lineal en K, tiene luego la representación:

4. La matriz de rigidez elemental:

i i K

1 en el nodo

0 en los nodos ,i

i

j k

i

( )w x

( ) ( ) ( ) ( ) ( ) ( ) ( )i i j j k kw x w N x w N x w N x

( ), ( ), ( ) base de fcs lineales en i j kx x x K

, , ,

, ,

sim. ,

i i i j i k

K j j j k

k k

a a a

a a

a

A


41

Ejemplo

( , )i i i ix y x y

, , /i i i ( , ) 1 1 1

( , ) 0 1 0

( , ) 0 1 0

i i i i i i

i j j j j i

i k k k k i

x y x y

x y x y

x y x y

Resolviendo:1

12 det 1 (área triángulo)2

1 ;2 2

j k k ji ii

j jj k k j

k ki i

x y x y x y

x yy y x x

x y

Luego:

2 2,i i

i i

i ii i i i i iK K K

i i

x x

y y

a d d d

x x x

,i j i j i j i jKa d x


42

Ejemplo

2 2

2 2

2 2sim.

i i i j i j i k i k

K j j j k j k

k k

A

Matriz de rigidez elemental:

Ilustraremos el proceso de ensamble para el caso siguiente:

Elemento 1) i=1; j=2; k=31 1 11 2 3

1 1 11 2 3

1 1 11 2 3

0 0 0

1 10

1 10

h h

h h

Sustituyendo:

11

1 2

3

2 1 11

1 1 02

1 0 1

A ξ


43

EjemploPor similaridad:

43

3 1

6

2 1 11

1 1 02

1 0 1

A ξ5

55 7

1

2 1 11

1 1 02

1 0 1

A ξ

Elemento 2) i=3; j=6; k=1 2 2 23 6 1

2 2 23 6 1

2 2 23 6 1

0 0 0

1 10

1 10

h h

h h

32

2 6

1

2 1 11

1 1 02

1 0 1

A ξ1

44 4

5

2 1 11

1 1 02

1 0 1

A ξ2

66 1

7

2 1 11

1 1 02

1 0 1

A ξ

Luego:


44

Ensamble primer elemento

1

2

3

4

5

6

7

2 1 1

1 1 0

1 0 11

2

A ξ

11

1 2

3

2 1 11

1 1 02

1 0 1

A ξ


45

Ensamble segundo elemento

1

2

3

4

5

6

7

1 1 0

1 2 1

0 1

2 1 1

1 1 0

1 0 11

2

1

A ξ

32

2 6

1

2 1 11

1 1 02

1 0 1

A ξ


46

Ensamble tercer elemento

1

2

3

4

5

6

7

2 1 1

1 1 0

1 0

1 1 0

1

11

2

1 1 0

1 2 1

0 1 1

2 1

0 1 1

A ξ

43

3 1

6

2 1 11

1 1 02

1 0 1

A ξ


47

Ensamble cuarto elemento

1

2

3

4

5

6

7

1 1 0

1 2 1

0 1

2 1 1

1 1 0

1 0 1

2 1 11 1 0

1 2 11

1

1 1 0

1 0 1

0 1

2

1

A ξ

14

4 4

5

2 1 11

1 1 02

1 0 1

A ξ


48

Ensamble quinto elemento

1

2

3

4

5

6

7

2 1 1

1 1 0

1

21 1 0

1 0 1

1 1 0

1 2 1

2 1

0 1 1

1 1 0

1

1 1

1

0

1 0

1 0

1

21

11

2 1

0 11

A ξ

55

5 7

1

2 1 11

1 1 02

1 0 1

A ξ


49

Ensamble sexto elemento

1

2

3

4

5

6

7

1 1 0

1 2 1

0 1 1

1 1 0

1 2 1

0 1 1

21 12 1 1

1 1 0

1 0

1 1 0

1 2 1

0

0

1 2 1

0 1 1

1 1

1 1 0

1

1

0

11

21

1

A ξ

26

6 1

7

2 1 11

1 1 02

1 0 1

A ξ


50

Sistema de ecuaciones global1 1

2 2

3 3

4 4

5 5

6 6

7 7

4 1 1 1 1 0 0

1 1.5 0 0 0 0 0.5

1 0 1.5 0 0 0.5 0

1 0 0 1.5 0 0.5 0

1 0 0 0 1.5 0 0.5

0 0 0.5 0.5 0 0.5 0

0 0.5 0 0 0.5 0 0.5

b

b

b

b

b

b

b

A ξ

Notar que la ecuación 1 está completa (no hay más elementos que contribuyan allí) :

1 2 3 4 5 14 b

2

1 2 3 4 5 61 6

hb f f f f f f

Puede mostrarse, de manera similar, que:

con: 1i i

if f d x

La ecuación es idéntica a la que se obtiene por diferencias finitas. El término a derecha cambia.


51

1. En formulaciones variacionales para solución de BVP, trabajaremos con espacios V más grandes que los vistos hasta ahora.

2. Sucesión de Cauchy: sucesión que verifica

La sucesión de Cauchy es convergente si cuando

3. Espacio de Hilbert: Es un espacio lineal V, equipado de producto escalar y norma asociada, completo (o sea que las sucesiones de Cauchy convergen respecto de la norma del espacio).

Espacios de Hilbert

0, / si ,i jN v v i j N

1 2, , , Vnv v v

0iv v i


52

Espacios de Hilbert: ejemplos 1D

• Espacio L2(I) de funciones de cuadrado integrable en I=(a,b)

dotado del producto interno:

y la norma:

La desigualdad de Cauchy:

Ejemplo:

22 I

L (I) : está definida en I y v v v dx

2 2L (I) L (I)( , )v w v w

1

122

2

2

L (I) I( , )v v dx v v

I

( , )v w vwdx

2

1( ) , I (0,1), L (I) si <

2v x x x v


53

Espacios de Hilbert: ejemplos 1D (cont.)

• Espacio H1(I) en I=(a,b):

dotado del producto interno:

y la norma:

• Espacio en I=(a,b):

• Dado el problema modelo:

luego:

• Notas: – es más grande que el espacio V de funciones lineales a trozos

que veníamos usando.– En MEF, la norma del error es simplemente la norma en H1(I).

1 10H (I) : H (I) y (a) (b) 0v v v v

1 12 2

1 1

22

H (I) H (I)I( , )v v v dx v v

1H (I) I( , )v w vw v w dx

10H (I)

12H (I) : , L (I)v v v

1 10 0 Hallar H (I)(V) / , , H (I)u u v f v v

'' ( ), I

(a)(

0)

(b)

u f x xD

u u

10H (I)


54

Espacios de Hilbert en d, d=2,3

1 12 2

1 1

2

H ( ) H ( )( , )v v v v dx v v

21

2: L ( ), L ( ), 1, , H ( )i

vv v i d

x

22 : está definida en y L ( ) v v v dx

1H ( )( , )v w vw v w dx

( , )v w vwdx

122v v dx

•

– Producto escalar:

– Norma:

•

– Producto escalar:

– Norma:

• 110 :H H( ) ( ), 0v v v


55

Espacios de Hilbert en d, d=2,3 (cont)

• Dado el problema modelo (ec. de Poisson + CB homogéneas)

luego:

o, equivalentemente:

• Nota: (V) se dice la formulación débil de (D), y la solución de (V) es la solución débil de (D). Ésta no es necesariamente una solución clásica de (D). Para que lo sea, u debe ser suficientemente regular de modo que u esté definida en el sentido clásico.

en

0 (

obre )

s

u f

uD

1 10 0

1 1 Hallar H ( ) / a( , ) ( , ) a( , ) ( , ) H ( )

2( )

2u u u f u v v f v vM

F( )u ( )F v

1 10 0 Hallar H ( ) H ( )( ) / u u v d f vV x v dx

a( , )u v ( , )f v


56

Interpretación geométrica del MEF

Consideremos el BVP de difusión-reacción con CB homogéneas:

luego:

en

0 sob)

re (D

u u f

u

1 10 0 Hallar H ( ) / H (( )) u u v dx uvV dx fv dx v

,u v ( , )f v

, 0 Vh hu u v v

Hallar V / ,( ) ( , ) Vh h h hh u u v f v vV

10H ( )Sea Vh un subespacio de dimensión finita de , ej. el espacio de funciones

lineales a trozos. El MEF aplicado al problema de difusión-reacción da:

10V H ( )hv Tomando en (V), (V)(Vh) resulta:

La solución MEF uh es la proyección con resp a , de la solución exacta u sobre Vh, o sea que uh es el elemento de Vh más próximo u con resp a ||.|| , i.e.:

uh

Vh

u

10Producto interno en H ( )

1 10 0H ( ) H ( )

Vh hu u u v v

10H ( )


57

CB naturales y esenciales

Consideremos el BVP de difusión-reacción con CB tipo Neumann:

luego:

o:

con

en

sob(

r)

e

u u f

ug

n

D

1H ( )u v uv dx fv dx gv dx v

,u v ( , )f v ( , )g v

1F( ) a( , ) ( , ) ( , )

2u u u f u g u

1 1 Hallar H ( ) / F( ) F( ) H( () ) u u vM v

1 Hallar H )) /( (V u

•La CB Neumann, que no tiene que ser impuesta explícitamente sobre u, se denomina CB natural.

•La CB Dirichlet u=u0 sobre , que debe ser satisfecha explícitamente por u, se conoce como CB esencial.


58

Problema continuo en forma abstracta

Objetivos:

• Dar un tratamiento unificado a muchos problemas de la Mecánica y la Física, a fin de no repetir el mismo argumento en distintos casos concretos.

• Entender la estructura básica del MEF.

Hipótesis: sea V un espacio de Hilbert con producto escalar (.,.)V y norma ||.||V, a(.,.) es una forma bilineal en V×V, y L(.) una forma lineal en V, tales que

1. a(.,.) es simétrica, i.e.,

2. a(.,.) es continua, i.e.,

3. a(.,.) es V-elíptica o coerciva, i.e.,

4. L(.) es continua, i.e.,

a( , ) a( , ) , Vu v v u u v

V V0 / a( , ) , Vu v u v u v

2

V0 / a( , ) Vv v v v

V0 / L( ) Vv v v


59

Formulación abstracta del MEF para problemas elípticos (cont)

V( ) F( )= min F

1 Hallar V / , con ( ) F( )= a( , ) ( )

2vu vM u v v vv L

Hallar V /( ) a( ,, L( ) V)V u v vu v

Teorema: (M) y (V) son equivalentes, i.e., u satisface (M) si y sólo si u satisface (V). Además, uV, y ésta verifica

V/ .u


60

Formulación abstracta del MEF para problemas elípticos (cont)

Demo.: sea vV y arbitrarios. Luego, v+uV, así que

Siendo

luego

Como g tiene un mínimo en 0, debe cumplirse

F( ) F( ) , u u v g tiene un mínimo e n =0

2

2

1g( ) a( , ) L( )

2

1a( , ) a( , ) a( , ) a( , ) L( ) L( )

2 2 2 2

1a( , ) L( ) a( , ) L( ) a( , )

2 2

u v u v u v

u u u v v u v v u v

u u u u v v v v

g (0) a( , ) ) )( (L 0u v v V

g(0) g( )

(QED)

g ( ) a( , ) L( ) a( , )u v v v v


61

Estimación del error de discretización

• Sea Vh un subespacio de V, con dimensión finita M, y sea {1, 2,…, M} una base para Vh, de modo que toda vVh puede representarse

• Usando Vh, obtenemos los problemas discretos análogos a (M) y (V):

– Como

– Como

– Forma matricial de (Vh): con

M

=1

= ,i i ii

v

Hallar V / F( ) F(( ) ), Vh hh h hM u u v v

V a( , ) L( ), 1,2, ,Mj h h j ju j

Hallar V ( / a( , )) L( ), Vh h hh hu u v v vV

M M

=1 1

V = , a( , ) L( ), 1,2, ,Mh h h i i i i j i ji i

u u j

Aξ b a( , ), L( ).ij i j j jA b


62

Estimación del error de discretización (cont)

• Dado que a(.,.) es simétrica:

la matriz A es simétrica

• Dado que a(.,.) es V-elíptica:

la matriz A es positiva definida

la matriz A es no singular

M M M M2

V1 1 1 1

a( , ) a( , ) a( , ) 0 si 0 ( )i i j j i i j ji j i j

v v v v

η Aη η 0

a( , ) a( , )i j j i ij jiA A

h! ! V solución de ( )h hu V ξ / Aξ = b


63

Estimación del error de discretización (cont)

• Teorema: Sea uV la solución de (V) y uhVhV. Entonces:

• Demo.:

1.

2.

3.

4.

5. Operando:

V VVh hu u u v v

solución de ( ) a( , ) L( ), V Vhu V u w w w solución de ( ) a( , ) L( ), Vh h h hu V u w w w

Restando: a( , ) 0, Vh hu u w w Sea = , con V arbitrario.h hw u v v

2

V

V V

a( , )

a( , ) a( , )

a( , )

a( , )

h h h

h h h

h h

h

h

u u u u u u

u u u u u u w

u u u u w

u u u v

u u u v

(QED)


64

Norma de energía

Considerando

1. a(.,.) es continua, i.e.,

2. a(.,.) es V-elíptica o coerciva, i.e.,

podemos definir una nueva norma llamada norma de energía:

Esta norma es equiv a ||.||V, i.e., constantes positivas , tal que

El producto escalar asociado a ||.||a es

La ec de error resulta

de donde

Medida en la norma de energía, uh es la mejor aproximación a u.

V V0 / a( , ) , Vu v u v u v

2

V0 / a( , ) Vv v v v

12

aa( , ) , Vv v v v

V a Vc C Vv v v v

c , C

a( , ) a( , ).u v u v

a( , ) 0 Vh hu u v v

a aVh hu u u v v


65

Ejemplos concretos

Ejemplo 1: sea

En este caso, es la forma débil del problema de Neumann

Se verifica

• a(.,.) es una forma bilineal simétrica

• a(.,.) es V-elíptica con =1

• a(.,.) es continua con =1

• L(.) es continua con

• Luego:

1 2V=H ( ), .

2

a( , )=

L( )= , L ( )

v w vw v w dx

v fv dx f

1

2

H ( )a( , )=v v v

en , 0 sob( r) e u

D u u fn

1a( , )=L( ), H ( )v w v v

1 12 2

1 1H ( ) H ( )a( , ) a( , ) a(w, )v w v v w v w

2L ( )f

2 2L ( ) L ( )L( )v fv dx f v

12H ( ) L ( )

u f

a( , )=a( , )v w w v

1 11

H ( ) H ( )H ( )hu u u v v


66

Ejemplos concretos (cont)

Ejemplo 2: sea

En este caso, es la forma débil del problema

Se verifica



• a(.,.) es V-elíptica con =1/2

– Demo.:


• Luego:

10 I I

V=H (I), I (0,1), a( , )= , L( )= .v w v w dx v fv dx

2 22

I I I

1

2v dx v dx v dx

en I, ( ) ( )( ) 0u f u a bD u

10a( , )=L( ), H (I)v w v v

1 12 2L (I) L (I) H (I) H (I)

a( , )v w v w v w

2L (I)f

2 2L (I) L (I)IL( )v fv dx f v

12H (I) L (I)

u f

a( , )=a( , )v w w v

1 110H (I) H (I)

2 H (I)hu u u v v

12

0

1 121

0 0 0 01 1 1 1 1

2 2 22 2

0 0 0 0 0

H (I) ( )= ( ) ( )

1

2

x x

v v x v y dy v x v dy v dy v dy

v dx v dx v dx v dx v dx


67


Ejemplo 3: sea

En este caso, es la forma débil del problema de Poisson

Se verifica



•

a(.,.) es V-elíptica con =1/(C+1), C tal que


• Luego:

1 20V=H ( ), , a( , )= , L( )= .v w v wdx v fv dx

2 Cv dx v v dx

en ,( 0 s) obre u uD f

10a( , )=L( ), H ( )v w v v

1

2 2

H ( )a( , )v v v v v v dx

2L ( )f

2 2L ( ) L ( )L( )v fv dx f v

12H ( ) L ( )

C 1u f

a( , )=a( , )v w w v

1 1

10H ( ) H ( )

C 1 H ( )hu u u v v

1 12 2L ( ) L ( ) H ( ) H ( )

a( , )v w v v v w


68


Ejemplo 4:

• Definimos los espacios

con norma

• La forma débil del problema (D) consiste en hallar tal que

• a(.,.) es una forma bilineal y simétrica


• a(.,.) es V-elíptica con =1/3


• Luego:

4

4 en I=(0,1),( (0) (0) (1) (1) )

d uf u uD u u

dx

12

2

2 22

H (I) I=v v v v dx

22H (I)= : , , L (I)v v v v

2 20H (I)= : H (I), (0) (0) (1) (1) 0v v v v v v

20I I

= H (I)

a( , ) L( )

u v dx fv dx v

u v v

20V=H (I)u

22H (I) L (I)

3u f 2 220H (I) H (I)

3 H (I)hu u u v v

2 22 2L (I) L (I) H (I) H (I)

a( , )v w v w v w

1

2 2 22

H (I) Ia( , )v v v v v v dx

2L (I)f


69


Ejemplo 5: Consideremos el problema bi-armónico:


con norma

• La forma débil del problema (D) se obtiene haciendo

• Se puede demostrar que a(.,.) es una forma bilineal, simétrica, continua y V-elíptica, así como L(.) es una forma lineal continua.

0 0

2

=L( ) =a( , )

( )

v u v

vfv dx u v dx u v dx u v ds u v dx u ds

n nV

2 2 22

22 2H ( )= : , , , , , L ( )

v v v v vv v

x y x y x y

2 20H ( )= : H ( ), 0 sobre

vv v v

n

2 en , 0 sob e ( ) ru

u f uDn

2

2 2 222 2 2 22

2 2H ( )=

v v v v vv v dx

x y x y x y


70


Ejemplo 6: Consideremos el problema estacionario de convección-difusión:

Supongamos ||/ moderado.


• L(.) es una forma lineal continua.

• a(.,.) es una forma bilineal, continua y V-elíptica, pero no simétrica.

• Teorema: si L(.) es una forma lineal continua, y a(.,.) es una forma bilineal, continua y V-elíptica, pero no simétrica, se puede demostrar que hay solución única a (V), y está acotada. Sin embargo, en este caso no existe problema de minimización asociado a (V).

0

=L( ) =a( , )

( )

v u v

ufv dx u u v dx u v u v dx v ds

nV

β β

2 en , , ,( ) 0 sobre u u fD u β β


71


Ejemplo 7: Consideremos el problema estacionario de conducción de calor en 3

Definimos el espacio


• L(.) es una forma lineal continua si f,gL2().• a(.,.) es una forma bilineal, simétrica, continua y V-elíptica si

11V : H ( ), 0 sobre v v v

2

2

a( , ) L( )

( )

u v v

fv dx u v dx u v dx u v dx u v dx g v ds

u v dx fv dx g v ds

V

k k k n k

k

1

2

en , Ecuación del calor

=0 en CB Dirichlet

en CB Neuman

( )

n

iju f k

u

g

D

u

k

k n

, / ( )ijc C c k C x x


72

Algunos espacios de elementos finitos

• Sea el dominio acotado representado

por la “triangulación” (o malla) de elementos finitos Th={K}.

– En 1D, el elemento K es un intervalo.

– En 2D, los elementos más comunes

son triángulos o cuadriláteros.

– En 3D, tetraedros o hexaedros.

• Los espacios Vh más comunes en elementos finitos consisten en funciones polinómicas por tramos definidas sobre la malla Th.

• La definición de un espacio de elementos finitos Vh requiere especificar

– La malla Th del dominio

– La naturaleza de las funciones vVh sobre cada elemento K (ej., lineal, cuadrática, cúbica, etc.)

– Los parámetros usados para definir dichas funciones.

K


73

Requisitos de regularidad

• BVP de 2º orden • BVP de 4º orden

• Si los espacios consisten de funciones polinómicas, resulta

2V H ( )h

1V H ( )h

2 1 1V H ( ) V C ( ) : , C ( ), 1, ,dh hi

vv v i

x

1 0V H ( ) V C ( ) : es continua en h h v v


74

Algunos ejemplos de elementos finitos en 2D

• Sea el dominio 2 con frontera poligonal .

• Sea Th={K} una triangulación de en triángulos K.


– P1(K) es el espacio de funciones lineales en K:

luego {1,x,y} es una base en P1(K) y dim P1(K)=3.

– P2(K) es el espacio de funciones cuadráticas en K:

luego {1,x,y,x2,xy, y2} es una base en P2(K) y dim P2(K)=6.

– En general:

2 200 10 01 20 11 02( , ) , ( , ) K, ijv x y a a x a y a x a xy a y x y a

P (K) : es un polinomio de grado en Kr v v r

00 10 01( , ) , ( , ) K, ijv x y a a x a y x y a

0

P (K) : ( , ) , ( , ) K,i jr ij ij

i j r

v v x y a x y x y a

( 1)( 2)dim P (K)

2r

r r

K

x x1º

x y2º


75

Elemento finito triangular lineal

• Sea el espacio de funciones lineales a trozos:

• Los parámetros necesarios para describir las funciones vVh se denominan grados de libertad (gdl) globales, y se eligen coincidentes con los valores de v en los nodos de la triangulación Th.

• Si KTh es un triángulo lineal de vértices (xi,yi), i=1,2,3, los gdl elementales son los valores de v en los dichos vértices.

01K

V : C ( ) y P (K), K Th hv v v


76

Elemento finito triangular lineal (cont)

1 1 1 2 1 3 1 1

2 2 1 2 2 3 2 2

3 3 1 2 3 3 3 3

ec (3.

( , )

( , ) 7

( , )

)

v x y c c x c y

v x y c c x c y

v x y c c x c y

1! P (K) / ( , ) , 1, 2,3i i iv v x y i

1 1 2 3P (K) ( , ) , ( , ) K, iv v x y c c x c y x y c

• Teorema: Sea KTh un triángulo de vértices (xi,yi), i=1,2,3. Una función vP1(K) está determinada de manera única por los gdl elementales, i.e., dados los valores i,

•Demo. 2: Notar que dim P1(K) = # gdl, i.e., # incógn = # ecs.

•Luego, det B0 implica que si vP1(K) y v(xi,yi)=0 para i=1,2,3, entonces

debe ser vº0.

•Esto se puede probar sin necesidad de calcular det B, tarea que se complica a

medida que se usan polinomios de mayor orden.

1 1

2 2

3 3

1

!solución para dados si det B= 1 0

1i

x y

x y

x y

det B= 2 area K 0 a b

1

3

2

K

a

b

(QED)

•Demo.:

Evaluando en los vértices

B es no singular !solución


77

Determinación de las funciones de base para el triángulo lineal

Toda función vP1(K) puede representarse

Las funciones de base i (º coord de área del triángulo K) verifican

lo que da lugar al sist de ecs

de donde

Análogamente:

3

1

( , ) ( , ) ( , ), ( , ) K.i i ii

v x y x y v x y x y

1 si ( , ) / ( , ) , 1, 2,3.

0 si i i i i i j j ij

i jx y a b x c y x y i j

i j

1 1 1 1 1 1 1 1

1 2 2 1 1 2 1 2

1 3 3 1 1 3 1 3

( , ) 1

( , ) 0

( , ) 0

x y a b x c y

x y a b x c y

x y a b x c y

1 1

2 2

3 3 2 3 3 21

1 1

2 2

3 3

1

0

0

1 2 area K

1

1

x y

x y

x y x y x ya

x y

x y

x y

1

2

3 3 21

1 1

2 2

3 3

1 1

1 0

1 0

1 2 area K

1

1

y

y

y y yb

x y

x y

x y

1

2

3 2 31

1 1

2 2

3 3

1 1

1 0

1 0

1 2 area K

1

1

x

x

x x xc

x y

x y

x y

1

31

2

1

1

3

1

2

2

1

3

1

2

3

3 1 1 32 2 area K

x y x ya

1 32 2 area K

y yb

3 1

2 2 area K

x xc

1 2 2 1

3 2 area K

x y x ya

2 13 2 area K

y yb

1 2

3 2 area K

x xc


78

Continuidad entre elementos triangulares lineales

• Dado

• Adoptando como gdl los nodos de la malla Th, podemos definir alternativamente

• Demo.: Para probar que (3.11)º(3.4), es necesario probar que la función vVh definida de acuerdo a (3.11) es continua no sólo en los nodos sino también a través de las fronteras interelementales.

– Sean K1 y K2 dos triángulos en Th que comparten el lado S y los nodos 1 y 2.

– Sea la restricción de v a Ki.

– v1=v2 en los nodos 1 y 2, y v1,v2 lineales.

v1=v2 a lo largo de todo el lado S

v es continua a través de S

1KV : P (K), K T , y es continua en los nodos ec (3.11)h hv v v

01K

V : C ( ) y P (K), K ec (3.3)Th hv v v

2

1

v2

K1

K2

v1

S

i1K

P (K )i iv v

0C ( )v (QED)


79

Elemento finito triangular cuadrático

• Sea el espacio de funciones cuadráticas a trozos:

• Sea KTh un triángulo de vértices xi=(xi,yi), i=1,2,3,

y sea xij=(xi+ xj)/2 el punto medio del lado ij, i<j, i,j=1,2,3.

• Teorema: toda función vP2(K) está únicamente determinada por los gdl

02K


12

3

2

1

13

23K

( ), 1, 2,3,

( ), , , 1, 2,3.i

ij

v i

v i j i j

x

x


80

Elemento finito triangular cuadrático (cont)

• Demo.: como dim P2(K)= #gdl =6, es suficiente probar que si vP2(K) y v(xi)=v(xij)=0 (con i<j, i,j=1,2,3), entonces debe ser vº0.

1. A lo largo del lado 23, v varía cuadráticamente y v=0 en x2, x3 y x23. Luego, v=0 x , y v puede escribirse

2. Ídem a lo largo del lado 13, luego

3. Evaluando en x12

23

1 2 0 0 1 2 1( ) ( ) ( ) , K, cte, , funciones de base en P (K).v w w x x x x

1 1 1 1 1 1( ) ( ) ( ), K, P (K), función de base en P (K).v w w x x x x

12 1 12 2 12 0 0 0

1 1( ) ( ) ( ) 0 0 0

2 2v w w w v x x x (QED)


81

Elemento finito triangular cuadrático (cont)

Toda función vP2(K) puede expresarse

Con las funciones de base en P2(K) dadas por

Es fácil verificar que i, i=1,…,6, conforman una base en P2(K) y además i(xj)=ij.

1 1 1

2 2 2

3 3 3

2 1

2 1

2 1

4 12º

3

2

1

6º13

5º23K

3 3 6

1 , 1, 1

2 1 ( ) 4 ( ) ( )i i i i j ij i ii i j i j i

v v v v

x x x

4 1 2

5 2 3

6 1 3

4

4

4

funciones asociadas a los nodos en el medio de los lados

funciones asociadas a los nodos en vértices


82

Continuidad entre elementos triangulares cuadráticos

• Dado

o, alternativamente,

• Demo.: Para probar que (a)º(b), es necesario probar que la función vVh definida de acuerdo a (b) es continua no sólo en los nodos sino también a través de las fronteras interelementales.

– Sean K1 y K2 dos triángulos en Th que comparten el lado S, con nodos 1, 2 y 4.


– v1=v2 en los nodos 1, 2 y 4, y v1,v2 cuadráticas.



2KV : P (K), K T , y es continua en los nodo ec s (b) h hv v v

02K

V : C ( ) y P (K), K ec (aT )h hv v v

i2K

P (K )i iv v

0C ( )v (QED) v2

v12

1

4

S


83

Elemento finito triangular cúbico

• Sea el espacio de funciones cúbicas a trozos:

• Sea KTh un triángulo de vértices xi=(xi,yi), i=1,2,3, y


123

( ), 1, 2,3

( ), , 1, 2,3,

( )

i

iij

v i

v i j i j

v

x

x

x

123 1 2 3

2 1, , 1,2,3,

3 31

3

iij i j i j i j

x x x

x x x x

03K


221

3

2

1223

K

112

332123

113331


84

Elemento finito triangular cúbico (cont)

• Demo.: como dim P3(K)= #gdl =10, es suficiente probar que si vP3(K) y v(xi)=v(xiij)=v(x123)=0 (con i,j=1,2,3, ij), entonces debe ser vº0.

1. v tiene variación cúbica a lo largo de los lados 12, 23 y 13, y v=0 en cuatro puntos de cada lado, luego v=0 en los tres lados y puede escribirse

2. Evaluando en x123

1 2 3 0 0 1 2 3 1( ) ( ) ( ) ( ) , K, cte, , , bases en P (K).v w w x x x x x

123 1 123 2 123 3 123 0 0 0

1 1 1( ) ( ) ( ) ( ) 0 0 0

3 3 3v w w w v x x x x

(QED)


85

Funciones de base para elemento finito triangular cúbico

1 1 1 1

13 1 3 2

2 4 1 2 1

93 1

2

10 1 2 327 5 221º

3

2

16 223º

K

4 112º

7 332º

10 123º

9 113º8 331º


86

Continuidad entre elementos triangulares cúbicos

• Dado



– Sean K1 y K2 dos triángulos en Th que comparten el lado S, de nodos extremos 1 y 2, y nodos intermedios 4 y 5.


– v1=v2 en los nodos 1, 2, 4 y 5, y v1,v2 cúbicas.



3KV : P (K), K T , y es continua en los nod s (b)o ec h hv v v

03K


i3K

P (K )i iv v

0C ( )v (QED)v2

v12

1

5

S

4


87

Elemento finito triangular cúbico con gdl en derivadas

• Sea el espacio de funciones cúbicas a trozos:

• Sea KTh un triángulo de vértices xi=(xi,yi),

i=1,2,3, y centro de gravedad x123.


3KV : P (K), K Th hv v

123

( ), 1, 2,3

( ), ( ), 1, 2,3

( )

i

i i

v i

v vi

x y

v

x

x x

x

3

2

1K 123


88

Elemento finito triangular cúbico con gdl en derivadas (cont)

• Demo.: como dim P3(K)= #gdl =10, es suficiente probar que si vP3(K) y

entonces debe ser vº0.

1. la función v tiene variación cúbica a lo largo del lado 12, y se anula así como su derivada en dos puntos del lado, luego v=0 en todos los puntos del lado 12.

2. Razonando análogamente con los lados 23 y 13, llegamos a

3. Evaluando en x123


123 1 123 2 123 3 123 0 0 0

1 1 1( ) ( ) ( ) ( ) 0 0 0

3 3 3v w w w v x x x x

123( ) ( ) ( ) ( ) 0, 1,2,3i i i

v vv v i

x y

x x x x


(QED)


89

Continuidad entre EF triangulares cúbicos con gdl en derivadas

• Dado



–Sean K1 y K2 dos triángulos en Th que comparten el lado S, de nodos 1 y 2.

–Sea la restricción de v a Ki.

–v1=v2 y (derivadas en la dirección s a lo largo de S) en los nodos

1 y 2, y v1,v2 cúbicas v1=v2 a lo largo de todo el lado S


• Nota: no se logra continuidad C1 por cuanto la función de base asociada a x123 no llega con

pendiente nula a los lados.

3KV : P (K), K T , y , , continuas en los n ec (bd )o osh h

v vv v v

x y

03K


i3K

P (K )i iv v

0C ( )v (QED)

1 2v v

s s


90

Elemento finito triangular C1-continuo

• Consideremos un espacio de EF , lo que requiere usar polinomios de grado 5 por triángulo, i.e.

• Sea KTh un triángulo de vértices xi=(xi,yi), i=1,2,3,

y sea xij=(xi+ xj)/2 el punto medio del lado ij, i<j, i,j=1,2,3.


1V C ( )h

15K


2 2 2

2 2

( ), 1, 2,3

( ), ( ), 1, 2,3

( ), ( ), ( ), 1, 2,3

( ), , 1, 2,3,

i

i i

i i i

ij

v i

v vi

x y

v v vi

x x y x

vi j i j

n

x

x x

x x x

x

12

3

2

1

13

23

n

K


91

Elemento finito triangular C1-continuo (cont)

• Demo.: como dim P5(K)= #gdl =21, es suficiente probar que si todos los gdl son nulos, entonces debe ser vº0.

1. v es un polinomio de grado 5 a lo largo del lado 23 y

2. es un polinomio de grado 4 a lo largo del lado 23 y

3. Aplicando idéntico razonamiento sobre los lados 12 y 13, llegamos a

23( ) ( ) ( ) 0, 2,3 0, 23i i

v v v vi

n n s n n

x x x x

12

3

2

1

13

23

s

K

n

2

2( ) ( ) ( ) 0, 2,3 0, 23i i i

v vv i v

s s

x x x x

v

n

2

1 3 3 30, 0, 23 ( ) ( ) ( ), K, P (K)v

v v p pn

x x x x x

2 2 2

0 1 2 3 0( ) ( ) ( ) ( ) , K, ctev w w x x x x x

5 0P K 0 ( ) 0, Kv w v º x x (QED)


92

Continuidad entre EF triangulares C1-continuos

• Sean K1 y K2 dos triángulos en Th que comparten el lado S, de extremos 2, 3.

• Sea la restricción de v a Ki., y w=v1v2 sobre S.

• Luego:

15K

5K

2 2 2

2 2

V ( ) : P (K), K T

: P (K), K T ,

, , , , , continuas en los nodos vértices

continua en los nodos medios de los lados

h h

h

v C v

v v

v v v v vv

x y x x y y

v

n

º

1C ( )v

(QED)

i5K

P (K )i iv v

23( ) ( ) ( ) 0, 2,3 0, Si i

w w w wi

n n s n n

x x x x

0, Sw

s

x

2

2( ) ( ) ( ) 0, 2,3 0, Si i i

w ww i w

s s

x x x x


93

Elemento finito tetraédrico lineal

• Sea la unión de un conjunto Th={K} de tetraedros no superpuestos K tales que ningún vértice de algún tetraedro se ubique sobre el lado de otro tetraedro.

• Adoptamos los siguientes espacios polinómicos por trozos de EF

• Para r=1, toda función vP1(K) está únicamente

determinada por sus valores en los vértices de K.

• En este caso, el espacio de EF es1

2

3

4

rV : P (K), K T , i.e., , ( , , ) K,i j mh h ijm ijm

i j m r

v v v a x y z x y z a

01V : C ( ) y P (K), K Th hv v v

K


94

Elemento finito rectangular bilineal

• Sea la unión de un conjunto Th={K} de rectángulos no superpuestos K tales que ningún vértice de algún rectángulo se ubique sobre el lado de otro rectángulo.

• Definimos el espacio

• Toda función vQ1(K) está únicamente determinada por sus valores en los vértices del rectángulo K.

• Se puede demostrar fácilmente que existe continuidad interelementos de v.

• Luego, el espacio de EF es

1 00 10 01 11Q (K) : es bilineal en K, i.e., ( , ) , ( , ) K, ijv v v x y a a x a y a xy x y a

1 2

34

01V : C ( ) y Q (K), K Th hv v v

K


95

Elemento finito rectangular bicuadrático

• Definimos el espacio de funciones bicuadráticas en K

• Toda función vQ2(K) está únicamente determinada por sus valores en los vértices, en el medio de los lados y en el centro del rectángulo K.

• Se puede demostrar fácilmente que existe continuidad interelementos de v.

• Luego, el espacio de EF es

2

2, 0

Q (K) : ( , ) , ( , ) K, i jij ij

i j

v v x y a x y x y a

02V : C ( ) y Q (K), K Th hv v v

K

1 2

34

5

6

7

8 9


96

Resumen

• Definimos un elemento finito como la terna {K, PK, }, donde

– K es un objeto geométrico.

– PK es un espacio lineal de dimensión finita de funciones definidas en K.

es un conjunto de gdl que determinan de manera única toda función vK.

• Por ej., para el EF triangular lineal {K, PK, }, resulta

– K es un triángulo.

– PK=P1(K).

es el conjunto de valores de v en los vértices de K.


97

Tipos de elementos finitos más comunes

Valor de la función

Valores de las derivadas primeras

Valores de las derivadas segundas

Valor de la derivada normal al lado

3 P (K) C1

06 P (K) C2

0

10 P (K) C3

0

10 P (K) C3

0

GeometriaGrados de libertad


# g

dl

# g

dl

Esp

aci

o d

efu

ncio

nes

Esp

aci

o d

efu

nci

on

es

Co

ntin

uid

ad

de

le

spa

cio

ME

F

Con

tinu

ida

d d

el

esp

aci

o M

EF

21 P (K) C5

1

18 P (K) C5'

1


98

Tipos de elementos finitos más comunes (cont)

Valor de la función

Valores de las derivadas primeras

Valores de las derivadas segundas

Valor de la derivada normal al lado

3 P (K) C2

0

10 P (K) C2

0

4 P (K) C3

1

2 P (K) C1

0

4 P (K) C1

0

4 Q (K) C1

0

9 Q (K) C2

0

16 Q (K) C3

0



# g

dl

# g

dl

Esp

aci

o d

efu

ncio

nes

Esp

acio

de

fun

cio

ne

s

Co

ntin

uid

ad

de

le

spa

cio

ME

F

Con

tinu

ida

d d

el

esp

aci

o M

EF


99

Soporte de diferentes funciones de base

Soporte de funciones de base asociadas a nodos de vértice

Soporte de funciones de base asociadas a nodos sobre los lados

Soporte de funciones de base asociadas a nodos en el centro


100

Estimación de error en problemas elípticos

• Para un típico problema elíptico de la forma

donde se verifica

1. a(.,.) es una forma bilineal simétrica, continua y V-elíptica.

2. L(.) es una forma lineal continua

resulta

• Si elegimos v=huV como un interpolante de u y estimamos el error de interpolación ||uhu||V, obtendremos una estimación del error ||uuh||V del MEF.

• Elegimos hu tal que sus gdl coincidan con los de u en Vh, así el problema de determinar ||uhu||V se reduce a determinar uhu individualmente sobre cada elemento finito KTh.

Hallar V / a( , )=L( ) Vu u v v v

V VVhu u u v v


101

Interpolación con funciones lineales a trozos en 2D

• Sea

• Para el triángulo KTh, definimos

• /hK da una idea de la calidad del elemento (cuanto mayor, mejor)

• Designemos Th a una familia de mallas {Th} caracterizadas por el parámetro , y una constante +, independiente de h, tal que

Esta condición implica que los triángulos KTh no pueden ser arbitrariamente finos. La constante es una medida del ángulo más pequeño para cualquier KTh.

1 1

KV=H ( ), V = V: P (K), K Th hv v

K

K

: diámetro de K = lado más largo de K

: diámetro del círculo inscrito en K

h

KK Tmax

h

h h

hK

K

hK

K

K

K

K Thh


102

Interpolación con funciones lineales a trozos en 2D (cont)

• Sean Ni, i=1,2,…,M, los nodos de Th. Dado , definimos el interpolante huVh por

i.e, hu es la función lineal a trozos que coincide con u en los nodos xi de Th.

• Empecemos por estimar el error uhu en cada triángulo K.

( ) ( ) 1, 2, ,i ihv v i M x x

0Cu

K

u hu


103


• Teorema: sea KTh un triángulo de vértices xi, i=1,2,3. Dado vC0(K), sea el interpolante vP1(K) definido por

Luego:

1.

2.

donde

• Nota: la magnitud de los errores de interpolación en la función y sus derivadas primeras dependen del valor de las derivadas segundas, que es una medida de cuán curva es la superficie descrita por la función.

( ) ( ), 1, 2,3.i iv v i x x

2KL (K) L (K)2

2 max Dv v h v

2K

L (K) L (K)1 2K

max D 6 max Dh

v v v

1 2 1 2 1 2D , , , .

vv

x y

L (K) Kmax ( )v v

xx


104


• Demo.: como vP1(K), usando las funciones de base i podemos escribir

• Usando una expansión de Taylor, en el punto y = x+ tenemos

• Tomando y=xi:

• Como ||xxi||hK xK, i=1,2,3, el resto Ri(x) resulta acotado por

2

1

( ) ( ) ( ) ( , )j jj j

vv v x y R

x

y x x x y

3 3

1 1

( ) ( ) ( ) ( ) ( )i ii i

i i

v v v

º x x x x x

22

, 1

1( , ) 0 1

2 i i j ji j i j

vR x y x y

x x

x y x Δ

2

1

( ) ( ) ( ) ( , )i i ij j

j j

vv v x x R

x

x x x x x

( )iR x( )ip x

2K L (K)2

( ) 2 max D K, 1,2,3iR h v i

x x


105


• Luego:

• Por otro lado, veremos más adelante que:

3

1

3 3

1

2K L (K)

1

2

( ) ( ) ( ) ( )

( )( ) ( )

2 max D K

( ) ( ) ( ) max ( )

i ii

i i i i ii

i i

v

v

v R

R R Rv

h v

x

x x x x

xx x x

x

x x

3 3

1 1

1, 0.i i ii i

p

3 3 3 3

1 1 1 1

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )ii i i i i i

i i i i

v v v p R

x x x x x x x x x

(QED 1)


106


• Se calcula la derivada de v:

• Por otro lado, veremos más adelante que:

3

1

3 3 3

K1 1 1K

2K

L (K)2K

1( ) max ( ) ( )

6 max D K

ii

ij j j

i ii i i

i i ij jj j

v v

x

v vR

x x x

x

h

R R

v

Rx x

x

x

x

x x

3 3

1 1

0, .i ii

i ij j j

vp

x x x

3 3 3 3

1 1 1 1

( )ii i i ii i

i i i ij j j j j

vv v v p R

x x x x x

x

(QED 2)

K1/

Alberto Cardona, Víctor Fachinotti Cimec-Intec (UNL/Conicet), Santa Fe, Argentina Introducción al...

Documents

Transcript of Alberto Cardona, Víctor Fachinotti Cimec-Intec (UNL/Conicet), Santa Fe, Argentina Introducción al...