Resolucion mediante el MEF de la ecuacion de
Laplace
Motivacion:
Problema estacionario de conduccion del calor.
Balance de la energa termica:
Tasa de cambio Calor generado Flujo de calor
de energa termica = por fuentes + entrante en
en un cuerpo en un cuerpo el cuerpo
Energa termica:
Z
c(x)(x)(x, t)dV
c(x) : capacidad calorfica
(x) : densidad
(x, t) : temperatura
Calor generado por fuentes:
Z
f (x, t)dV
f (x, t) : densidad de calor generado
Flujo de calor entrante:
Figura 1: Dominio de estudio y normal saliente
1
Z
q ndS
q : calor
: contorno de
n : normal saliente de
Balance:
Z Z Z
d
c(x)(x)(x, t)dV = f (x, t)dV q ndS
dt
Por el teorema de la divergencia:
Z Z
q ndS = qdV
V
Para el caso estacionario:
Z
d
c(x)(x)(x, t)dV = 0
dt
Por tanto, Z Z
f (x, t)dV q(x, t)dV = 0
Z
[f (x, t) q(x, t)] dV = 0
Esto se debe cumplirse para cualquier por pequeno que sea. Por tanto
q = f en
Ley de Fourier (ecuacion constitutiva)
q = k
k : coeficiente de conductividad termica
Por tanto,
(k) = f en
k = f en
donde, para R2 ,
2 2
= = , (Laplaciano)
x21 x22
Entonces, se obtiene la siguiente expresion:
2
k = f en
Si k = 1 y f = 0,
= 0 en (ecuacion de Laplace)
Problema de valores de contorno
= 0 en (ecuacion de gobierno)
= en D (condiciones de Dirichlet)
= () n = q en N (condiciones de Neumann)
n
Forma debil
Multiplicando ambos miembros de la ecuacion de Laplace e integrando sobre
el dominio , se obtiene
Z
wdV = 0 w
Por otro lado, se tiene que
(w) = w + w
= w + w
Entonces,
w = (w) w
Reemplazando lo anterior, se obtiene
Z Z Z
wdV = (w)dV w dV = 0 w
Z Z
w dV n (w)dS = 0 w
Se consideran funciones de peso que cumplen w = 0 en D y, por tanto,
Z Z
w dV w(n )dS = 0 w
N
Z Z
w dV wq dS = 0 w
N
donde en la ultima expresion se han considerado las condiciones de contorno de
Neumann. A esta expresion se la conoce como forma debil de la ecuacion de
Laplace.
3
De la forma debil a la forma fuerte
Se tiene que
Z Z Z
w dV = (w)dV wdV w
Z Z
= n (w)dS wdV w
Z Z
= (n )wdS wdV w
N
Por tanto, teniendo en cuenta que n = n ,
Z Z Z
wdV + wdS wq dS = 0 w
N n N
y Z Z
wdV + q wdS = 0
w
N n
Considerese el subespacio W0 = {w : w | = 0} del espacio de funciones ad-
misibles. Entonces,
Z
wdV = 0 w W0
Considerese una funcion > 0 en y = 0 en . Adoptese w = W0 .
Entonces Z
()2 dV = 0
= 0 en
Con lo que se obtiene
Z
q wdS = 0 w
N n
q = 0 en N
n
Aproximacion de la forma debil
A partir de
Z Z
w dV wq dS = 0 w
N
Se construye el problema aproximado
Z Z
wh h dV wh q dS = 0 wh
N
4
con
h (x) = 1 N1 (x) + . . . + n Nn (x)
wh (x) = w1 N1 (x) + . . . + wn Nn (x)
Se puede ver facilmente que
X
h (x) = 1 N1 (x) + . . . + n Nn (x) = i Ni (x)
i
X
wh (x) = w1 N1 (x) + . . . + wn Nn (x) = wj Nj (x)
j
Entonces,
!
Z Z X X
h h
w dV = wj Nj (x) i Ni (x) dV
j i
Defnanse
N1 N1
x1 x2 w1
1
.. .. .. ..
T
B =
. .
;
w= . ; = .
Nn Nn
wn
n
x1 x2
Entonces,
wh = Bw
h = B
y
(wh ) (h ) = wT B T B
Definiendo Z
K= B T BdV
se obtiene que Z
wh h dV = wT K
5
Espacios de elementos finitos
La malla de elementos finitos
Dado un dominio poligonal se crea una particion de elementos {1 , 2 , ..., e , ..., E },
donde
E : numero de elementos Esto es, e f = 0 para e 6= f
[E
e = , =
e=1
Se considera e poligonal.
Se pide que cada lado de e sea parte del contorno o lado de algun otro
elemento. Esto deja fuera a configuraciones de elementos como la mostrada en
la figura (2).
Figura 2: Configuracion de elementos no permitida
Puntos nodales
Estan ubicadas por lo menos en los vertices del elemento, como se indica en
la figura (3).
Figura 3: Puntos nodales de un elemento triangular de tres nodos
Pero pueden ubicarse en otros puntos.
Existen siempre n nodos: 1, 2, ..., n cuyas posiciones son x1 , x2 , ..., xn
El conjunto de nodos y elementos forman la malla de elementos finitos.
6
Funciones base Ni
1. Son continuas y acotadas
Ni C()
2. Existen n funciones Ni , una por cada nodo
N1 , N2 , ..., Nn
y
Ni (x) = 0 x e
si
xi 6 e
Esto es, Ni es cero en los elementos que no estan conectados al nodo i
3. Ni (xj ) = ij
4. Ni (x) es un polinomio en e si xi e
Observacion:
(e) (e)
Ni (xej ) = ij en e con Ni (x) = NT L(e,i) (x)
Aproximacion de la forma debil con elementos finitos
Z Z Z E Z
X
wh h dV = wh h dV +. . .+ wh h dV = wh h dV
1 E e=1 e
donde Z
wh h dV = wT K e
e
con
Z
Ke = = B T BdV
e
Considerando una aproximacion por elementos finitos, se tiene que, para e ,
(e) (e)
N1 Nn
0 ... x1
...
x1
... 0
B=
(e) (e)
N1 Nn
0 ... ... ... 0
x2 x2
donde n(e) : denota el numero de nodos de e
Defnase Z
Ke = = BeT Be dV
(e)
7
con (e) (e)
N1 Nn
x ...
1 x1
Be =
N (e) Nn
(e)
1
...
x2 x2
Entonces Z
wh h dV = weT K e e
con
(e) (e)
w1 1
.. ..
we = . ; e = .
w(e)
(e)
n n
Por otro lado,
(e)
Z n
X Z
h
w q dS = wh q dS
e e
Donde e = e
Z Z
wh q dS = wT N q dS
e e
con
0
..
.
(e)
N1 (x)
..
N= .
(e)
Nn (x)
..
.
0
Defnanse, Z
fs(e) = N (e) q dS
e
y
(e)
N1
..
N (e) = .
N (e)
n
Entonces, Z
wh q dS = weT fs(e) en e
e
8
Globalmente,
Z Z
h h
w dV wh q dS = 0 w
Implica
wT (Ku fs ) = 0 w
(e)
Observacion. Para encontrar K y fs hay que ensamblar K (e) y fs
Elementos triangulares de tres nodos
Figura 4: Elemento triangular de tres nodos
Coordenadas triangulares
Figura 5: Coordenadas triangulares
Observese la figura (5). A partir de ella se puede definir la coordenada tri-
angular asociada al vertice i como
A i = 0 para los nodos j 6= i
i =
A i = 1 para el nodo i
9
Figura 6: Grafica de la coordenada triangular
Se tiene, entonces, un sistema de coordenadas (1 , 2 , 3 ) En este sistema,
las coordenadas de los nodos son
Nodo 1: (1, 0, 0)
Nodo 2: (0, 1, 0)
Nodo 3: (0, 0, 1)
Por otro lado, las doordenadas del centroide son 31 , 31 , 13
Las coordenadas 1 , 2 , 3 forman una particion de la unidad:
1 + 2 + 3 = 1
Interpolacion lineal en R2
Si f (x, y) vara linealmente, en e , entonces f (x, y) = a0 + a1 x + a2 y en e .
Por otro lado, existen tres funciones lineales tales que
1 = N1 (x, y)
2 = N2 (x, y)
3 = N3 (x, y)
Defnase, ademas,
f1 = f (x1 , y1 )
f2 = f (x2 , y2 )
f3 = f (x3 , y3 )
De lo anterior se sigue que
f (x, y) = f1 N1 (x, y) + f2 N2 (x, y) + f3 N3 (x, y)
o, de manera equivalente,
N1 (x, y)
f (x, y) = f1 f2 f3 N2 (x, y)
N3 (x, y)
10
Transformacion de coordenadas
Se tiene que
1 + 2 + 3 = 1
Ademas,
1 x1 + 2 x2 + 3 x3 = x
1 y1 + 2 y2 + 3 y3 = y
Por tanto,
1 1 1 1 1
x = x1 x2 x3 2
y y1 y2 y3 3
Se puede demostrar que
1 1 1
2A = det x1 x2 x3
y1 y2 y3
El cambio de coordenadas inverso es
1 x y x3 y2 y2 y3 x3 x2 1
1 2 3
2 = x3 y1 x1 y3 y3 y1 x1 x3 x
2A
3 x1 y2 x2 y1 y1 y2 x2 x1 y
De aqu,
i Ni (x, y)
= = yj yk j 6= k, i 6= k, i 6= k
x x
i Ni (x, y)
= = xk xj j 6= k, i 6= k, i 6= k
y y
Se puede entonces plantear las siguiente aproximacion lineal para la incognita
dentro de un elemento:
(e) (e) (e) (e) (e) (e)
h (x) = u1 N1 (x) + u2 N2 (x) + u3 N3 (x) en e
(e) (e)
= u1 1 + u2 2 + u3 1 en e
(e)
En este caso, la matriz B tiene la siguiente forma:
(e) (e) (e)
N1 N2 N3
B (e) = x x x = 1 y2 y3 y3 y1 y1 y2
N (e) N (e) N (e) 2A x3 x2 x1 x3 x2 x1
1 2 3
y y y
Por tanto, Z T T
K (e) = B (e) B (e) dV = B (e) B (e) Ae
e
donde Ae es el area del elemento. Se obtiene finalmente que
Z T
wh h dV = w(e) K
e
11