Apuntes de la Clase de Elementos Finitos
Elaborado por: Felipe Fernández, PhD
Part I
Introducción a la Teorı́a de Elementos
Finitos
1 ¿Qué es el método de elementos finitos ?
Es un método numérico que resuelve un modelo fı́sico que usualmente se encuentra definido
por ecuaciones diferenciales parciales con sus respectivas condiciones iniciales y de frontera.
El método de elementos finitos se caracteriza por una formulación variacional, una es-
trategia de discretización, algoritmos de solución, y procedimientos de procesamiento pos-
terior, cf. Figura 1.
Problema estructural
Modelización
Formulación variacional
Discretización
Solución
Postprocesamiento
error? e Interpretación
óptimo?
Fin
Figure 1: Proceso de solución de un problema estructural usando el método de elementos
finitos.
2 Conceptos preliminares
El objetivo de esta sección es entender los terminos matemáticos que usamos para describir
los modelos fı́sicos en ecuaciones diferenciales. Por ejemplo, la ecuación fundamental de este
1
curso es la del equilibrio que se la puede escribir ası́:
σij,j (x) + bi (x) = 0 x ∈ Ω. (1)
Describiremos a continuación conceptos que nos ayuden a entender esta y otras ecuaciones.
2.1 Conjuntos (sets)
Colección de objetos con ciertas caracteristicas. Ejemplos:
• Conjunto de números reales, R := {x : −∞ < x < +∞}
• Conjunto de números positivos, R+ := {x ∈ R : x > 0}
• ρ-vencindario de x, Nρ (x) := {y ∈ Rn : |y − x| ≤ ρ}
• dominio Ω, conjunto de puntos x ∈ Rn que son parte de la estructura
• dominio Ωh , conjunto de puntos x ∈ Rn que son parte de la estructura discretizada
(Ωh ⊆ Ω)
2.2 Vectores
Vector x ∈ Rn es un arreglo de n números, dónde a n se le conoce como la dimensión del
vector. El vector columna de dimensión n podemos representarlo de la siguiente manera
x1
x2
x=
... (2)
xn
donde xi son las componentes del vector y el subindice puede tomar valores 1 < i < n. Otra
forma de representar un vector es con sus vectores base êi :
x = x1 ê1 + x2 ê2 + . . . + xn ên . (3)
Los vectores base son
1 0 0
0 1
, · · · ên = 0 .
ê1 =
. . . , ê2 = . . .
... (4)
0 0 1
La adición de dos vectores a, b ∈ Rn está definido como la suma componente por componente
a + b = (a1 + b1 )ê1 + (a2 + b2 )ê2 + · + (an + bn )ên . Ejemplos de vectores son la posición,
desplazamiento, velocidad, y las fuerzas.
Propiedades de los vectores:
• a+b=b+a
• (a + b) + c = a + (b + c)
• Existe 0 tal que a + 0 = a
• Existe un vector negativo −a tal que −a + a = 0
• α(βa) = (αβ)a
• (α + β)a = αa + βa
• 1a = a
2
2.3 Función
Regla que toma elementos del conjunto X (dominio) y produce un único elemento del
conjunto Y (codominio), f : X → Y, y = f (x).
Función lineal cumple f (αa + βb) = α f (a) + β f (b).
Los argumentos de una función son las variables que se los identifica usando los
paréntesis. Los argumentos se diferencias con los parámetros que son fijos o constantes.
Por ejemplo, supongamos que tenemos la ecuación x2 − d = 0
• r : R → R, r(x) = x2 − d (d es parámetro fijo)
• r̃ : R × R → R, r̃(x, d) = x2 − d (d es variable)
• r̂ : R×R → R, r̂(x̂(d), d) = (x̂(d))2 −d (x̂ es variable que es función de la variable
d)
2.4 Derivadas
Considerando una función escalar cuyo argumento es un vector f : x ∈ Rm → R, la derivada
y gradiente se definen como
D f (x) = [∂f /∂x1 , ∂f /∂x2 , · · · , ∂f /∂xm ] , (5)
∂f /∂x1
∂f /∂x2
∇ f (x) =
···
∂f /∂xm
= ∂f /∂x1 ê1 + ∂f /∂x2 ê2 + · · · + ∂f /∂xm êm . (6)
El operador gradiente toma un escalar y lo convierte en vector y por ello se lo puede definir
como ∇ = ∂/∂x1 ê1 + ∂/∂x2 , ê2 + · · · + ∂/∂xm êm .
2.5 Producto punto o escalar o interior
Dados a, b ∈ Rn , el producto escalar se lo define y representa
n
X
a·b= ai bi = a⊤ b = ai bi . (7)
i=1
Ejemplo del producto punto es la definción de trabajo como fuerza producto punto desplaza-
miento w = f · d.
Propiedades del producto punto:
• a·b=b·a
• (αa + βb) · c = αa · c + βb · c
• a·a≥0
• Si a · a = 0, entonces a = 0
√
• La norma o magnitud de un vector es definida |a| := a·a
• a y b son perpendiculares, si a · b = 0
• Ángulo θ entre dos vectores a y b está dado por cos θ = a·b
|a| |b|
• êi · êj = δij , la delta de Kronecker es definida por δij = 1 si i = j, y δij = 0 si i ̸= j
3
2.6 Espacio de vectores
Un conjunto V es un espacio vectorial real, cuyos elementos son vectores, si está equipado
con dos operaciones: suma y una multiplicación escalar.
2.7 Producto cruz o vectorial
El resultado de un producto cruz entre dos vectores a y b es otro vector que es perpendicular
a ambos. En tres dimensiones, el producto vectorial de las bases son
ê1 × ê2 = +ê3 , ê2 × ê3 = +ê1 , ê3 × ê1 = +ê2 ,
ê2 × ê1 = −ê3 , ê3 × ê2 = −ê1 , ê1 × ê3 = −ê2 , (8)
ê1 × ê1 = 0 , ê2 × ê2 = 0 , ê3 × ê3 = 0 .
Con estas definiciones podemos calcular cualquier producto tal que
a × b = (a1 ê1 + a2 ê2 + a3 ê3 ) × (b1 ê1 + b2 ê2 + b3 ê3 )
= a1 ê1 × (b1 ê1 + b2 ê2 + b3 ê3 ) +
a2 ê2 × (b1 ê1 + b2 ê2 + b3 ê3 ) +
a3 ê3 × (b1 ê1 + b2 ê2 + b3 ê3 )
:0
= a1 b1ê1×ê1 + a1 b2 ê1 × ê2 + a1 b3 ê1 × ê3 +
:0
a2 b1 ê2 × ê1 + a2 b2 ×ê2 + a2 b3 ê2 × ê3 +
ê2
:0
a3 b1 ê3 × ê1 + a3 b2 ê3 × ê2 + a3 b3ê3×ê3
= a1 b2 ê3 − a1 b3 ê2 − a2 b1 ê3 + a2 b3 ê1 a3 b1 ê2 − a3 b2 ê1
= (a2 b3 − a3 b2 )ê1 + (a3 b1 − a1 b2 ê3 )ê2 + (a1 b2 − a2 b1 )ê3
ê1 ê2 ê3
= a1 a2 a3 . (9)
b1 b2 b3
La magnitud del producto cruz entre dos vectores esta dado por |a × b| = |a||b| sin(θ) donde
θ es el ángulo entre los dos vectores. Ejemplo de producto cruz es el torque pues esta dado
por el producto cruz entre la fuerza y la posición t = f × x.
2.8 Producto diádico o tensorial
Definimos el producto diádico como el producto matricial
⊤
a1 b1 a1 b1 a1 b2 · · · a1 bn
a 2 b2
a2 b1 a2 b2 · · · a2 bn
a ⊗ b = ab⊤ =
= . (10)
... ... ··· ··· ··· ···
an bn an b1 an b2 · · · an bn
Por ejemplo, en dos dimensiones tenemos el producto diadico de los vectores base
1 0 0 1
ê1 ⊗ ê1 = , ê1 ⊗ ê2 = ,
0 0 0 0
(11)
0 0 0 0
ê2 ⊗ ê1 = , ê2 ⊗ ê2 = .
1 0 0 1
Propiedades del producto tensorial
4
• (αa + βb) ⊗ c = α(a ⊗ c) + β(b ⊗ c)
• a ⊗ (αb + βc) = α(a ⊗ b) + β(a ⊗ c)
• (a ⊗ b) c = (c · b) a
• (a ⊗ b) (c ⊗ d) = (b · c) (a ⊗ d)
2.9 Tensores
A un tensor de segundo órden (o de rango 2) lo podemos entender como un vector de vectores
o una matriz; por ello podemos usar el producto diádico de las bases de la siguiente forma
para definirlo
A11 A12 · · · A1n
A21 A22 · · · A2n
A= ··· ··· ··· ···
An1 An2 · · · Ann
= A11 ê1 ⊗ ê1 + A12 ê1 ⊗ ê2 + · · · A1n ê1 ⊗ ên +
A21 ê2 ⊗ ê1 + A22 ê2 ⊗ ê2 + · · · A2n ê2 ⊗ ên +
···
An1 ên ⊗ ê1 + An2 ên ⊗ ê2 + · · · Ann ên ⊗ ên
= Aij êi ⊗ êj . (12)
Ejemplos de tensores de segundo órden son: esfuerzo o estrés, y deformación. El tensor de
segundo órden en 1, 2 y 3 dimensiones seria
A11 A12 A13
A11 A12
A1d = A11 , A2d = , A3d = A21 A22 A23 . (13)
A21 A22
A31 A32 A33
Propiedades de los tensores de segundo órden
• A+B=B+A
• (A + B) + C = A + (B + C)
• O+A=A
• (−A) + A = O
• α(βA) = (αβ)A
• (α + β)A = αA + βA
• α(A + B) = αA + αB
• 1A = A
• (A + B)a = Aa + Ba
• (αA)a = α(Aa)
• Oa = 0
5
• Ia = a
• α(AB) = (αA)B = A(αB)
• A(B + C) = AB + AC
• (A + B)C = AC + BC
• A(BC) = (AB)C
• AO = OA = O
• IA = AI = A
• A⊤
ij = Aji
• a · (Ab) = (A⊤ a) · b
• (A ⊗ B)C = (C · B)A
Asociar un tensor con una matriz es ventajoso computacionalmente. Considere los sigu-
ientes ejemplos
a = Ab
ai êi = (Aij êi ⊗ êj ) (bk êk )
= (Aij bk ) (êi ⊗ êj ) êk
= (Aij bk ) (êj · êk ) êi
= (Aij bk ) δjk êi
= (Aik bk ) êi (14)
A B = (Aij êi ⊗ êj ) (Akl êk ⊗ êl )
= (Aij Bkl ) (êi ⊗ êj ) (êk ⊗ êl )
= (Aij Akl ) (êj · êk ) (êi ⊗ êl )
= (Aij Akl ) δjk (êi ⊗ êl )
= Aik Akl (êi ⊗ êl ) (15)
2.10 Divergencia
El operador divergencia sobre un tensor de órden 2 en 3D está definido
∂A11 ∂A12 ∂A13
+ +
∂x211 ∂x2 ∂x3
div A(x) = ∂A ∂A22 ∂A23
∂x1 + ∂x2 + ∂x3
∂A31 ∂A32 ∂A33
∂x1 + ∂x2 + ∂x3
∂Aij
= êi
∂xj
= Aij,j êi (16)
También podemos calcular la divergencia de un vector en 3D
∂a1 ∂a2 ∂a3
div a(x) = + +
∂x1 ∂x2 ∂x3
∂ai
=
∂xi
= ai,i . (17)
6
Por ejemplo, teniendo un campo vectorial en dos dimensiones f (x) = x31 ê1 + x32 ê2 , la
divergencia está dada por div (f (x)) = 3x21 + 3x32 . En el origen div (f (0)) = 0, mientras que
un punto xa = [1 , 1]⊤ tenemos div (f (xa )) = 6.
Propiedades de la divergencia
• div (α a) = α div (a)
• div (a) = ∇ · a
• div (A) = ∇ · (A⊤ )
• Si tensor es simétrico A = A⊤ , div (A) = ∇ · A
• div (f a) = f div (a) + (∇f ) · a
2.11 Teorema de la divergencia
Considerando que la integral de una derivada es la función evaluada en su frontera:
Z xb
∂f (x)
dx = f (xb ) − f (xa ) . (18)
xa ∂x
Podemos generalizar esto con el teorema de la divergencia
Z Z
div (f ) dv = f · n da , (19)
Ω Γ
donde Γ es la frontera superficial del volumen Ω, y n es el vector normal a la superficie de
esta frontera Γ. Similarmente, tenemos la siguiente expresión para un tensor
Z Z
div (A) dv = A n da . (20)
Ω Γ
El teorema de la divergencia lo usaremos muchas veces para pasar de una integral de
volumen hacia una de superficie en nuestro curso de elementos finitos. Pero existe ocasiones
que la integral de volumen es mas fácil de calcular, por ejemplo: Si tenemos el campo
5
vectorial
R f (x) = (3x1 + x3 )ê1 + (x22 − sin(x21 x3 ))ê2 + (x1 x3 + x2 ex1 )ê3 y deseamos calcular
Γ f · n da en la superficie de un cubo dado por la region 0 ≤ x1 ≤ 1, 0 ≤ x2 ≤ 3, 0 ≤ x3 ≤ 2,
7
para ello mejor usamos el teorema de la divergencia de la siguiente manera
Z Z
f · n da = div (f ) dv
Γ ZΩ
= (3 + 2x2 + x1 ) dv
Ω
Z 1 Z 3 Z 2
= (3 + 2x2 + x1 ) dx3 dx2 dx1
x1 =0 x2 =0 x3 =0
Z 1 Z 3
= (3 + 2x2 + x1 ) x3 |2x3 =0 dx2 dx1
x1 =0 x2 =0
Z 1 Z 3
= (6 + 4x2 + 2x1 ) dx2 dx1
x1 =0 x2 =0
Z 1
6x2 + 2x22 + 2x1 x2 |3x2 =0 dx1
=
x1 =0
Z 1
= (36 + 6x1 ) dx1
x1 =0
= 36x1 + 3x21 |1x1 =0
= 39 . (21)
2.12 Tensor de cuarto órden
Ası́ como un tensor de segundo órden actuando en un vector resulta en otro vector Ab = c,
podemos generalizar que un tensor de cuarto órden actuando en un tensor de segundo órden
resulta en otro tensor de segundo órden
AB = C, (22)
donde A es un tensor de cuarto órden. A este también se lo puede definir con sus bases
(similar a la Ecuación (12)):
A = A1111 ê1 ⊗ ê1 ⊗ ê1 ⊗ ê1 + A1112 ê1 ⊗ ê1 ⊗ ê1 ⊗ ê2 + · · ·
A111n ê1 ⊗ ê1 ⊗ ê1 ⊗ ên + · · ·
Annn1 ên ⊗ ên ⊗ ên ⊗ ê1 + Annn2 ên ⊗ ên ⊗ ên ⊗ ê2 + · · ·
Annnn ên ⊗ ên ⊗ ên ⊗ ên
= Aijkl êi ⊗ êj ⊗ êk ⊗ êl . (23)
Cuando tenemos elasticidad lineal, un ejemplo de tensor de cuarto órden es la rigidez o
tensor de elasticidad C. Este tensor actuando sobre las deformaciones obtenemos el tensor
de estrés, i.e., σ = C ϵ.
Propiedades de los tensores de cuarto órden
• A(αB + βC) = αAB + βAC
• (A + B)C = AC + BC
• (αA)B = α(AB)
• OA = O
• IA = A
8
• A+B=B+A
• A+O=A
• A + (−A) = O
• α(βA) = (αβ)A
• (α + β)A = αA + βA
• 1A = A
• (α + β)A = αA + βA
• α(A + B) = αA + αB
• (A + B)C = AC + BC
• OA = O
• IA = A
3 Problema Elástico
tp
Γt
Γ
Γu
x
ê2
ê1
0
Figure 2: Problema elástico con su dominio y condiciones de frontera.
Un problema elástico lineal queda definido por:
• 3 ecuaciones de equilibrio (cons sus condiciones de frontera)
• 6 ecuaciones constitutivas
• 6 ecuaciones de compatibilidad.
3.1 Ecuaciones de equilibrio
La ecuación de equilibrio está dada por:
∂σij (x)
+ bi (x) = 0 , x ∈ Ω, (24)
∂xj
donde σij son las componentes del estrés, y bi son las componentes de la fuerza de volumen.
El ejemplo mas común de fuerza de volumen de nuestro interés es la gravedad dado por
9
b = −ρ g ê2 donde ρ y g son la densidad y la aceleración de la gravedad respectivamente.
Las ecuaciones de equilibrio en realidad son tres:
∂σ11 (x) ∂σ12 (x) ∂σ13 (x)
+ + + b1 (x) = 0 , x∈Ω
∂x1 ∂x2 ∂x3
∂σ21 (x) ∂σ22 (x) ∂σ23 (x)
+ + + b2 (x) = 0
∂x1 ∂x2 ∂x3
∂σ31 (x) ∂σ32 (x) ∂σ33 (x)
+ + + b3 (x) = 0 . (25)
∂x1 ∂x2 ∂x3
3.2 Condiciones de frontera
Sobre la superificie Γu , nosotros prescribimos las condiciones de frontera de los desplaza-
mientos, también llamadas esenciales o Dirichlet:
u(x) = up para x ∈ Γu . (26)
En la superficie complementaria Γt , prescribimos las condiciones de frontera de las tracciones,
también llamadas naturales o Neumann:
σ(x) n(x) = σ p
para x ∈ Γt . (27)
S T
Que estas superficies sean complementarias significa que Γu Γt = Γ y Γu Γt = ∅.
3.3 Ecuaciones de compatibilidad
El tensor de pequeñas deformaciones se lo calcula usando la gradiente de los deplazamientos:
1
ϵ= ∇u + (∇u)⊤
2
1 ∂ui ∂uj
= + êi ⊗ êj
2 ∂xj ∂xi
∂u1 1 ∂u1 ∂u2 1 ∂u1 ∂u3
+ +
∂x1 2 ∂x2 ∂x1 2 ∂x3 ∂x1
= ··· ∂u 1 ∂u ∂u (28)
∂x2
2
2 ∂x3
2
+ ∂x2
3
∂u3
Sym ··· ∂x3
Por simetrı́a, tenemos solo 6 ecuaciones.
3.4 Ecuaciones constitutivas
Como vimos, el tensor de rigidéz o elasticidad C es de cuarto órden y transforma las defor-
maciones en estreses, cf. σ = Cϵ. Por las simetrias tanto del tensor estrés y deformación,
podemos simplificar y poner estas ecuaciones de manera mas compacta. Para esto vector-
izamos y usamos una notación reducida de los tensores de estrés y deformaciones, conocida
como la convención de Voigt:
[σ1 , σ2 , σ3 , σ4 , σ5 , σ6 ] = [σ11 , σ22 , σ33 , σ12 , σ23 , σ13 ], (29)
[ϵ1 , ϵ2 , ϵ3 , ϵ4 , ϵ5 , ϵ6 ] = [ϵ11 , ϵ22 , ϵ33 , 2ϵ12 , 2ϵ23 , 2ϵ13 ] . (30)
Note que en esta notación los coeficientes de los tensores de estrés y deformacion van de 1 al
6, por ser componentes independientes de cada entidad. Con esta notación nuestro tensor
10
de rigidez o elasticidad de manera general (material anisotrópico) se la puede representar
con matriz:
σ1 ϵ1 c11 c12 c13 c14 c15 c16 ϵ1
σ2 ϵ2 · · · c22 c23 c24 c25 c26 ϵ2
= [C] ϵ3 = · · · · · · c33 c34 c35 c36
σ3
ϵ3 .
ϵ4 · · · (31)
σ4
· · · · · · c44 c45 c46 ϵ4
σ5 ϵ5 · · · · · · · · · · · · c55 c56 ϵ5
σ6 ϵ6 Sym · · · ··· ··· · · · c66 ϵ6
Estas representan 6 ecuaciones. Note que los coeficientes reducidos corresponden a los
coeficientes del tensor de cuarto órden, por ejemplo c11 = C1111 y c34 = C3312 . Para un
modelo isotrópico elático lineal, la forma matricial del tensor de elasticidad depende sólo de
dos coeficientes. Estos coeficientes pueden ser las constantes de Lamé (λ y µ) o el módulo
de Young E y un coeficiente de Poisson ν:
λ + 2µ λ λ 0 0 0
··· λ + 2µ λ 0 0 0
··· ··· λ + 2µ 0 0 0
[C] = , (32)
··· ··· ··· µ 0 0
··· ··· ··· · · · µ 0
Sym ··· ··· ··· ··· µ
Eν E
λ= (1+ν)(1−2ν) , µ=G= 2(1+ν) . (33)
4 Tarea 1
⊤ ⊤ 1 2 4
1. Considere los siguientes vectores y matrices x = [1 , 2] , y = [5 , 3 , 8] , A = ,
6 2 1
2 3
B= . Con ellos calcule las siguientes expresiones
5 1
(a) x · x
(b) y ⊗ y
(c) x y⊤
(d) A y
(e) x · (A y)
(f) A⊤
(g) A⊤ x · y
(h) A B
(i) B A =
(j) A⊤ B
⊤
(k) A⊤ B⊤
2. De la pregunta 1, ¿cuál es la dimensión de cada vector?
3. De las matrices de la pregunta 1, ¿cuál(es) podrı́an ser considerados tensores? Y si
son tensores, ¿cuál es su dimensión y su rango?
4. Evalue la gradiente para las siguientes funciones, asuma que x ∈ R2
11
(a) f (x) = |x|
(b) f (x) = (4x1 + 9x2 )/(x1 x2 )
5. Demuestre que si un tensor A es simétrico, entonces div A = ∇ · A
6. Calcule la divergencia de los siguientes campos vectoriales:
(a) f (x) = x21 x2 ê1 − (x33 − 3x1 )ê2 + 4x22 ê3
x31 x22
(b) f (x) = (3x1 + 2x23 )ê1 + x3 ê2 − (x3 − 7x1 )ê3
7. Siendo el dominio un cuboide 0 R≤ x1 ≤ 1, 0 R≤ x2 ≤ 3, 0 ≤ x3 ≤ 2, calcule por
separado las siguientes integrales Ω div (f ) dv y Γ f ·n da para las funciones vectoriales
del Ejercicio 6a y compruebe que el teorema de la divergencia se cumple.
8. Proponga una función del desplazamiento sencilla cuyas pequeñas deformaciones re-
sulten en un tensor constante con deformaciones sólo en las componentes 11, 22, 23 y
32. ¿Cuál serı́a su estrés si el material es isotrópico (use las constantes de Lamé por
conveniencia)?
9. Calcule las tracciones superficiales
en las 6 caras
de un cuerpo cúbico si estuviera con
2 4 5
un estado de estrés dado por: σ = 4 6 3 mN2
5 3 9
10. ¿Es el siguiente campo de estrés en equilibrio estático? σ11 = cx1 x2 , σ12 = σ21 =
c(a2 − x22 ), σ22 = σ33 = σ23 = σ13 = 0. Si no, ¿Cuál deberia ser la fuerza del cuerpo?
12
