DESARROLLO DE UN SOFTWARE EN

ELEMENTOS FINITOS COMPATIBLE CON

GID PARA RESOLUCIÓN DE
PROBLEMAS ELÁSTICO LINEALES.

Verónica Gallego Otalvaro

Coordinador Santiago Correa
Grupo de Bioingeniería
EAFIT
10/12/2009
Contenido

INTRODUCCIÓN .................................................................................................3

MARCO TEÓRICO ..............................................................................................4

1. Elemento triangular lineal bidimensional ................................................................................. 6

2. Elemento triangular cuadrático bidimensional ..................................................................... 12

3. Elemento rectangular lineal bidimensional ............................................................................ 19

4. Elemento rectangular cuadrático .............................................................................................. 24

5. Elementos tridimensionales ....................................................................................................... 30

6. Tetraedro lineal............................................................................................................................... 31

7. Hexaedro lineal ............................................................................................................................... 36

8. Integración numérica .................................................................................................................... 40

9. Ensamble de la matriz de rigidez (K) ........................................................................................ 42

METODOLOGÍA ................................................................................................50

TRABAJO FUTURO..........................................................................................50

CONCLUSIONES ..............................................................................................51

REFERENCIAS .................................................................................................52

2
Introducción

El método de los elementos finitos es un método numérico utilizado frecuentemente

para obtener la solución aproximada de un sistema de ecuaciones diferenciales
parciales, permitiendo encontrar soluciones apropiadas para diversos capos de la
ingeniería, uno de los tópicos en el que se utiliza mas comúnmente es en los
problemas físicos que poseen geometrías muy complicadas lo que hace necesario el
uso de un computador para llevar a cabo la implementación del método.

El método es altamente utilizado ya que permite obtener una solución numérica

aproximada sobre el comportamiento de una estructura bien sea un sólido, un gas o
un fluido, sobre el que están definidas ciertas ecuaciones diferenciales en forma débil
o integral que se encargan de definir el comportamiento físico del problema
dividiéndolo en un número elevado de subdominios no intersectantes entre sí, esto se
denomina elementos finitos.

Dentro de cada elemento se encuentran una serie de puntos que permiten el

movimiento de cada elemento, estos puntos son denominados nodos y son
adyacentes entre ellos si pertenecen al mismo elemento, también un nodo sobre la
frontera puede pertenecer a varios elementos, el conjunto de nodos y las relaciones
que existen entre ellos constituyen una malla, esta malla es de gran importancia ya
que a partir de esta se realiza la discretización del dominio en elementos finitos, y a su
vez permite realizar los cálculos sobre cada uno de los nodos. La construcción de la
malla generalmente se realiza a partir de programas especializados, pero en este libro
se mostrara una forma de realizar la implementación en un lenguaje de programación
como lo es Matlab, para esto se hace necesario entender las relaciones de adyacencia
o conectividad que poseen las variables que no conocemos en cada uno de los nodos
que además nos permiten saber los grados de libertad que posee la pieza a estudiar.
El conjunto de relaciones entre el valor de una determinada variable y los nodos se
puede expresar como un sistema de ecuaciones donde el número de ecuaciones es
proporcional al número de nodos, este sistema también pede ser escrito como una
matriz llamada matriz de rigidez global.

Generalmente el método de elementos finitos se programa computacionalmente para

calcular los desplazamientos, las deformaciones y tensiones cuando se trata de un
problema de sólidos deformables o más generalmente un problema de mecánica de
medios continuos, además el método es fácilmente adaptable a problemas de
transmisión de calor, de mecánica de fluidos para calcular campos de velocidades y
presiones o de campo electromagnético. Dada la dificultad de encontrar la solución a
estos problemas de manera analítica, con frecuencia en la práctica los métodos
numéricos y, en particular, los elementos finitos, se convierten en una herramienta
alternativa para encontrar estas soluciones, ya que la convergencia de esta
metodologías puede decirse es buena.

3
Marco teórico

Para resolver un problema definido mediante ecuaciones diferenciales por medio de

un algoritmo de elementos hay que tener en cuenta las siguientes etapas.

1. El problema debe reformularse en forma variacional.

2. El dominio de variables independientes debe dividirse mediante una partición en

subdominios. Asociada a la partición anterior se construye un espacio vectorial
de dimensión finita, llamado espacio de elementos finitos. Siendo la solución
numérica aproximada obtenida por elementos finitos una combinación lineal en
dicho espacio vectorial.

3. Se obtiene la proyección del problema variacional original sobre el espacio de

elementos finitos obtenido de la partición. Esto da lugar aún sistema con un
número de ecuaciones finito, aunque en general con un número elevado de
ecuaciones incógnitas. El número de incógnitas será igual a la dimensión del
espacio vectorial de elementos finitos obtenido y, en general, cuanto mayor sea
dicha dimensión mejor será la aproximación numérica.

4. Cálculo numérico de la solución del sistema de ecuaciones.

Los pasos anteriores permiten construir un problema de cálculo diferencial en un

problema de álgebra lineal. Dicho problema en general se plantea sobre un espacio
vectorial de dimensión no-finita, pero que puede resolverse aproximadamente
encontrando una proyección sobre un subespacio de dimensión finita, y por tanto con
un número finito de ecuaciones. La discretización en elementos finitos ayuda a
construir un algoritmo de proyección sencillo, logrando además que la solución por el
método de elementos finitos sea generalmente exacta en un conjunto finito de puntos.

Formulación débil: La formulación débil de una ecuación diferencial permite convertir

un problema de cálculo diferencial formulado en término de ecuaciones diferenciales
en términos de un problema de álgebra lineal planteado sobre un espacio de Banach,
generalmente de dimensión no finita, pero que puede ser aproximado por un sistema
finito de ecuaciones algebraicas.

Dada una ecuación diferencial lineal de la forma:

ξ (u ) = f

Donde la solución es una cierta función definida sobre un dominio Ω ⊂ ℜ n , y se han

especificado un conjunto de condiciones de contorno adecuadas, puede suponerse
que la función buscada es un elemento de un espacio de funciones o espacio de
Banach V .

 u ∈V f ∈V
A(u ) = f , 
A :V → V ′ A ∈ ξ (V ,V ′)

4
Donde V ′ es el espacio dual de V , la forma variacional débil se obtiene buscando la
única solución u ∈ V tal que:

 a ( u , v ) = Au , v

a (u , v ) = f , v , ∀v ∈ V donde 
 f , v = ∫ fvdx
 Ω

Cuando el operador lineal es un operador elíptico, el problema se puede plantear

como un problema de minimización sobre el espacio de Banach.

Discretización del dominio: Dado un dominio Ω ⊂ ℜ n con una frontera continua en

el sentido de Lipschitz una partición en n "elementos finitos", es una colección de n
subdominios {Ω ( e ) }e =1 que satisface:
n

1. Ω ⊂ ℜ n
2. Cada Ω ( e ) es un conjunto compacto con una frontera Lipschitz-continua.
3. int(Ω ( i ) ) ∩ int(Ω ( j ) ) = 0, i≠ j

Usualmente por conveniencia práctica y sencillez de análisis, todos los "elementos

finitos" tienen la misma "forma", es decir, existe un dominio de referencia Ω̂ ⊂ ℜ n y una
colección de funciones biyectivas de la siguiente forma:

{F (e) ˆ → Ω (e)
F (e) : Ω }
En los análisis 2D el dominio de referencia Ω̂ se suele tomar como un triángulo
equilátero o un cuadrado, mientras que en los análisis 3D, el dominio de referencia
típicamente es un tetraedro o un hexaedro. Además sobre cada elemento se
considerarán algunos puntos especiales, llamados nodos y que generalmente incluirán
los vértices del elemento finito y se requerirá la condición adicional de que dos
elementos adyacentes tengan los compartan los nodos sobre el subconjunto
Ω ( i ) ∩ Ω ( j ) , es decir:

x ∈ Ω ( i ) ∩ Ω ( j ) ∧ ( x ∈ nod (Ω ( i ) )) → ( x ∈ nod (Ω ( j ) ))

Una vez definida la partición en elementos finitos, se define sobre cada elemento un
espacio funcional de dimensión finita, usualmente formado por polinomios de grado.
Este espacio funcional servirá para aproximar localmente la solución del problema
variacional. El problema variacional en su forma débil se plantea sobre un espacio de
dimensión no-finita, y por tanto la función buscada será una función de dicho espacio.
El problema en esa forma exacta es computacionalmente inabordable, así que en la
práctica se considerará un subespacio de dimensión finita V h del espacio vectorial
original V . Y en lugar de la solución exacta presentada anteriormente se automatizará
la proyección de la solución original sobre dicho el subespacio vectorial de dimensión
finita, es decir, se resolverá numéricamente el siguiente problema:

a (u h , v h ) = f , v h , ∀v h ∈ V h

5
Donde

u h = ∏ e ( u ) ∈ V h , es la solución aproximada.

∏e :V → V h V h ⊂ V Es el proyector ortogonal del espacio original sobre el

subespacio vectorial asociado a la discretiación.

Si la discretización es suficientemente fina y el espacio funcional finito sobre cada

elemento está bien escogido, la solución numérica obtenida aproximará
razonablemente bien la solución original. Eso implicará en general considerar un
número muy elevado de elementos finitos y por tanto un subespacio de proyección de
dimensión elevada. El error entre la solución exacta y la solución aproximada puede
acotarse gracias al lema de Ceá, que en esencia afirma que la solución exacta y la
solución aproximada satisfacen:

u − uh ≤ inf
h h
u − vh
v v ∈V v

Es decir, el error dependerá ante todo de lo bien que el subespacio vectorial asociado
a la discretización en elementos finitos V h aproxime el espacio vectorial original V .

1. Elemento triangular lineal bidimensional

1.1 Funciones de forma

Los elementos finitos pueden ser clasificados por grupos de acuerdo al orden
polinomial que posean, lineal, cuadrático, cúbico o de mayor orden.

Fig.1: elemento triangular lineal

La función de forma de un elemento triangular lineal bidimensional puede se escrita a

partir de un polinomio que posee elementos lineales y un termino constante como se
muestra en la ecuación (1.1)

φ = α1 + α 2 x + α 3 y (1.1)

6
El polinomio es lineal en x , en y , y además contiene unos coeficientes que hacen
referencia a cada uno de los nodos.

El polinomio de interpolación del elemento triangular lineal se define por la relación

(1.1), donde φ se utiliza para representar cualquier variable desconocida. En cada
nodo.

φ = φ1 para x = x1 y = y1
φ = φ2 para x = x 2 y = y2 (1.3)
φ = φ3 para x = x3 y = y3

Así,

φ = α 1 + α 2 x1 + α 3 y1
φ = α 1 + α 2 x2 + α 3 y 2 (1.4)
φ = α 1 + α 2 x3 + α 3 y 3

Ahora para encontrar los valores de α1 , α 2 y α 3 , tenemos:

1
α1 = [( x 2 y 3 − x3 y 2 )φ1 + ( x3 y1 − x1 y 3 )φ 2 + ( x1 y 2 − x 2 y1 )φ3 ]
2A
1
α2 = [( y 2 − y 3 )φ1 + ( y3 − y1 )φ 2 + ( y1 − y 2 )φ3 ] (1.5)
2A
1
α3 = [( x3 − x 2 )φ1 + ( x1 − x3 )φ 2 + ( x 2 − x1 )φ3 ]
2A

Donde el área, A , está definida por:

2 A = ( x1 y 2 − x 2 y1 ) + ( x3 y1 − x1 y 3 ) + ( x 2 y 3 − x3 y 2 ) (1.6)

El área del elemento está definida en términos de las coordenadas nodales.

Obteniendo las funciones de forma, donde se realiza una aproximación de φ con una
función que depende de las tres funciones nodales, como se presenta a continuación.

φ =  1φ1 +  2φ 2 +  3φ3 (1.7)

Las funciones de forma quedan definidas como

1
 1 ( x, y ) = [( x 2 y3 − x3 y 2 ) + ( y 2 − y3 ) x + ( x3 − x 2 ) y ] (1.8.1)
2A
1
 2 ( x, y ) = [( x3 y1 − x1 y3 ) + ( y3 − y1 ) x + ( x1 − x3 ) y ] (1.8.2)
2A
1
 3 ( x, y ) = [( x1 y 2 − x 2 y1 ) + ( y1 − y 2 ) x + ( x 2 − x1 ) y ] (1.8.3)
2A

7
Estas ecuaciones se ajustan a la relación general expresada por la ecuación (1.1).
Como ejemplo, el valor de 1 en el nodo 1 se puede obtener mediante la sustitución
de x = x1 , y = y1 en la ecuación (1.8.1) así obtenemos:

1
1 = [x2 y3 − x3 y 2 + y 2 x1 − y3 x1 + x3 y1 − x2 y1 ] (1.9)
2A

1 es igual a cero en los nodos 2, 3, y todos los puntos de una línea que pasan a
través de estos nodos. El valor entre paréntesis en la ecuación (1.9) también es igual a
2 A , por lo tanto tenemos

2A
1 = =1 con x = x1 , y = y1
2A
1 = 0 con x = x 2 , y = y2 y x = x3 , y = y3 (1.10)

Fig.2: Representación grafica de la ecuación (1.10) cuando 1 = 1

El gradiente de la variable φ está dado por la siguiente expresión

∂φ ∂ 1 ∂ ∂
= φ1 + 2 φ 2 + 3 φ 3
∂x ∂x ∂x ∂x
∂φ ∂ 1 ∂ 2 ∂ 3
= φ1 + φ2 + φ3 (1.11)
∂y ∂y ∂y ∂y

El valor de ∂1 ∂x , se puede determinar fácilmente de la ecuación (1.8.1)

∂ 1 y 2 + y 3
= (1.12)
∂x 2A

El resto de las derivadas de las funciones de forma se obtienen de una manera

sencilla y el valor de la derivada se convierte en

∂φ 1
= [( y 2 − y 3 )φ1 + ( y 3 − y1 )φ 2 + ( y1 − y 2 )φ3 ] (1.13)
∂x 2 A

8
Sin embargo es algo complicado encontrar los α i por medio de este sistema de
ecuaciones por lo cual se hace necesario escribir las funciones de forma en términos
de las coordenadas de área.

1.2 Coordenadas de área

A fin de simplificar el proceso de solución utilizando elementos triangulares, se

construye un sistema de coordenadas bidimensionales que permite facilitar la
evaluación de las integrales que se presentaran posteriormente

Se construyen las coordinadas de área definidas a partir de las variables, L1 , L2 y L3

A1
L1 =
A
A
L2 = 2 (1.17)
A
A
L3 = 3
A

Cuando A1 , A2 y A3 son áreas parcialmente definidas por la unión de un punto P en

el triángulo con los nodos en las esquinas (Fig.3), de manera tal que Ai es el área
situada frente al nodo i .

Fig.3: Sistema coordenado de área

Las líneas Li serán paralelas al lado del triángulo opuesto al nodo i . Teniendo en
cuenta que L1 + L2 + L3 = 1 , en cualquier punto dentro del triángulo. La Fig.4 muestra
la dirección cada vez mayor para cada coordenada. Es evidente que cada Li varía
entre 0 y el largo del lado opuesto el nodo i .

9
Para un triángulo lineal, podemos definir las funciones de forma como:

1 = L1
 2 = L2 (1.18)
 3 = L3

Fig.4: Coordenadas de área en un triangulo lineal

Por lo tanto, las coordenadas ( x, y ) están relacionadas con las coordenadas

( L1 , L2 , L3 ) a través de las siguientes expresiones

x = L1 x1 + L2 x2 + L3 x3
x = L1 y1 + L2 y 2 + L3 y3 (1.19)
1 = L1 + L2 + L3

Los cálculos de las derivadas de las funciones de forma requieren algunas

modificaciones ya que las coordenadas no son independientes, es decir,
L1 + L2 + L3 = 1 . Lo que se hace habitualmente es asumir L1 y L2 son
independientes. Así, utilizando la regla de la cadena tenemos

∂ i ∂ i ∂x ∂ i ∂y
= +
∂L1 ∂x ∂L1 ∂y ∂L1
∂ i ∂ i ∂x ∂ i ∂y
= + (1.21)
∂L2 ∂x ∂L2 ∂y ∂L2

La ecuación (1.21) puede ser reescrita como

 ∂ i   ∂x ∂y   ∂ i   ∂ i 
 ∂L   ∂L ∂L1     
 1= 1  ⋅  ∂x  = J  ∂x  (1.22)
 ∂ i   ∂x ∂y   ∂ i   ∂ i 
 ∂L   ∂L2 
∂L2   ∂y   ∂y 
 2

10
La matriz de 2 x 2, se llama la matriz Jacobiana (J). Además se pude ver que L3 es
una función de φ

∂φ ∂φ ∂L1 ∂φ ∂L2 ∂φ ∂L3

= + + (1.23)
∂L1 ∂L1 ∂L1 ∂L2 ∂L1 ∂L3 ∂L1

Sin embargo, L1 y L2 se supone que son independientes, y L3 = 1 − L1 − L2 por lo

tanto, se tiene que:

L3
= −1 (1.24)
L1

En consecuencia, las derivadas de φ respecto a las coordenadas naturales son las

siguientes:

∂φ ∂φ ∂φ
≡ −
∂L1 ∂L1 ∂L3
∂φ ∂φ ∂φ
≡ − (1.25)
∂L2 ∂L2 ∂L2

A continuación se presenta el código en Matlab par la construcción de la matriz B y el

área A de un elemento triangular de 3 nodos.
%-----------------------------------------------------------------------------
% Calculo Del Área
function A = calcular_A(coordenate_node1,coordenate_node2,coordenate_node3)

matrix=[ 1 coordenate_node1;...
1 coordenate_node2;...
1 coordenate_node3;];

A = 0.5*det(matrix);
end

%-----------------------------------------------------------------------------
% Calculo de B
function B = calcular_B(A,coordenate_node1,coordenate_node2,coordenate_node3)

b1 = coordenate_node2(2)-coordenate_node3(2);
b2 = coordenate_node3(2)-coordenate_node1(2);
b3 = coordenate_node1(2)-coordenate_node2(2);
c1 = coordenate_node3(1)-coordenate_node2(1);
c2 = coordenate_node1(1)-coordenate_node3(1);
c3 = coordenate_node2(1)-coordenate_node1(1);

B = (1/(2*A))* [b1 0 b2 0 b3 0 ;...

0 c1 0 c2 0 c3;...
c1 b1 c2 b2 c3 b3];
end
%-----------------------------------------------------------------------------

11
2. Elemento triangular cuadrático bidimensional

2.1 Funciones de forma

La función de forma de un elemento triangular cuadrático bidimensional puede se

escrita a partir de un polinomio que posee elementos lineales, de orden superior y un
termino constante. Sin importar que la forma sea igual a un elemento lineal
bidimensional, la interpolación polinomial de un elemento cuadrático es:

φ = α 1 + α 2 x + α 3 y + α 4 x 2 + α 5 xy + α 6 y 2 ( 2.1)

Los seis coeficientes indican que el elemento debe tener seis nodos.

Los elementos cúbicos o de orden superior también utilizan polinomios con términos
lineales y de orden superior, sin embargo, los límites del elemento deben ser paralelos
a los ejes de coordenadas.

Se puede realizar una aproximación de segundo grado sobre el elemento triangular.

El elemento triangular bidimensional tiene seis nodos, como se muestra en la Fig. 5.

Los α se determinan de la misma manera que para el elemento lineal. Definiendo φ

por cada nodo.

φ1 = α 1 + α 2 x1 + α 3 y1 + α 4 x12 + α 5 x1 y1 + α 6 y12
φ 2 = α 1 + α 2 x2 + α 3 y 2 + α 4 x22 + α 5 x2 y 2 + α 6 y 22 (2.2)


φ6 = α 1 + α 2 x6 + α 3 y 6 + α 4 x62 + α 5 x6 y 6 + α 6 y 62

Fig.5: Elemento triangular cuadrático

12
Como se menciono anteriormente esta metodología es un poco complicada para
encontrar los α por lo cual se escriben las ecuaciones de forma con respecto a las
coordenadas de área de la misma manera que se realizo para los triángulos.

Para el triángulo de segundo orden que se muestra en la Fig.5, las funciones de forma,
se escriben como funciones de L1 , L2 , y L3 :

1 = L1 (2 L1 − 1)
 2 = L2 (2 L2 − 1)
 3 = L3 (2 L3 − 1)
 4 = 4 L2 L3
 5 = 4 L1 L3
 6 = 4 L1 L2 (2.3)

Como el objetivo es encontrar la matriz de rigidez para este elemento es necesario

construir la matriz B y encontrar el jacobiano a partir de las funciones de forma,
además en este caso no es tan fácil resolver las integrales que se presentan para
encontrar la matriz K (matriz de rigidez), como en el caso del triangulo lineal por lo que
se hace necesario realizar una integración numérica.

Ejemplo 1: con el fin de ilustrar el uso de las coordenadas de área, se encontrara el

valor de ∂ 5 ∂y y ∂ 5 ∂x en el punto (1,1) en el triángulo que se muestra en figura 6

Fig.6: Ejemplo para calcular el valor de ∂ 5 ∂y y ∂ 5 ∂x

Los valores de x y y están definidos por la siguiente expresión.

x = L1 x1 + L2 x2 + L3 x3
y = L1 y1 + L2 y2 + L3 y3 (2.4)

Sustituyendo por los valores nodales obtenemos.

x = 3L2 + L3
y = L2 + 3L3 (2.5)

13
La derivada de la ecuación (1.25) utilizando la matriz Jacobiana, se convierten en:

∂x ∂x ∂x
≡ − = −1
∂L1 ∂L1 ∂L3
∂x ∂x ∂x
≡ − =2
∂L2 ∂L2 ∂L2
∂y ∂y ∂y
≡ − = −3
∂L1 ∂L1 ∂L3
∂y ∂y ∂y
≡ − = −2 ( 2. 6)
∂L2 ∂L2 ∂L2

El Jacobiano de la ecuación (1.22) es

 − 1 − 3
J =  (2.7)
 2 − 2

Como queremos encontrar la ∂ 5 ∂y y ∂ 5 ∂x un términos del nodo 5, reordenamos

la ecuación (1.22) en la siguiente forma

 ∂ i   ∂ i 
 ∂x  −1  
 ∂  = J  ∂∂x  (2.8)
 i  i
 ∂y   ∂x 

Ahora determinemos ∂ 5 ∂y y ∂ 5 ∂x para el nodo 5,

∂ 5 ∂ 5 ∂ 5
= − = 4 L3 − 4 L1
∂L1 ∂L1 ∂L3
∂ 5 ∂ 5 ∂ 5
= − = −4 L1 (2.9)
∂L2 ∂L2 ∂L2

Sustituyendo en la ecuación (2.9)

 ∂ i 
 ∂x  −1 
4 L3 − 4 L1 
 ∂  = J   (2.10)
 i  − 4 L1 
 ∂y 

Donde la inversa del Jacobiano esta dada por

1 − 2 3 
J −1 =  (2.11)
8 − 2 − 1

Reemplazando el Jacobiano en la ecuación (2.11)

14
  ∂ i 
 ∂x  1 − 2 3  4 L3 − 4 L1 
 ∂  =    
 i  8 − 2 − 1  − 4 L1 
 ∂y 
1  − 8 L3 − 4 L1 
=
8 − 8L3 + 12 L1 

 1 
 − L3 − 2 L1 
=  (2.12)
3 
− L3 + L1 
 2 

Estas dos ecuaciones son aplicables en todo el elemento. Para encontrar los valores
reales de las dos derivadas, determinamos L1 , L2 y L3 en el punto (1,1) mediante la
ecuación (4.32) y el hecho de que L1 + L2 + L3 = 1

1 = 0 L1 + 3L2 + L3
1 = 0 L1 + L2 + 3L3
1 = L1 + L2 + L3 (2.13)

Resolviendo el sistema de ecuaciones tenemos

1
L1 =
2
1
L2 =
4
1
L3 = (2.14)
4

Los valores finales de las dos derivadas son

∂ 5 1 1
= − L3 − L1 = −
∂x 2 2
∂ 5 3 1
= − L3 + L1 = (2.15)
∂y 2 2

A continuación se presenta el código en Matlab par la construcción de la matriz B de

un elemento triangular de 6 nodos.
%-----------------------------------------------------------------------------
%Calcular B
%Esta función entrega la matriz B y su determinante y lo que requiere como
%entradas son los puntos x, y y los puntos de integración que se explicaran
%mas adelante.
function [B detJ1] = calcular_B(x,y,L,p)

% Calculo de las derivadas de x e y respecto a los ejes coordenados

dxdl1 = x(1)-x(3);
dxdl2 = x(2)-x(3);

15
 dydl1 = y(1)-x(3);
dydl2 = y(2)-x(3);

% Calculo de las derivadas de las funciones de forma respecto a las

% coordenadas parametrizadas
dN1dl1 = 4*L(p,1)-1;
dN1dl2 = 0;
dN2dl1 = 0;
dN2dl2 = 4*L(p,2)-1;
dN3dl1 = 0;
dN3dl2 = 0;
dN4dl1 = 0;
dN4dl2 = 4*L(p,3);
dN5dl1 = 4*L(p,3);
dN5dl2 = 0;
dN6dl1 = 4*L(p,2);
dN6dl2 = 4*L(p,1);

% Almacenamiento del jacobiano

J1 = [dxdl1 dydl1; ...
dxdl2 dydl2];

detJ1 = det(J1);

% Almacenamiento de la matriz de derivadas respecto a coordenadas

%parametrizadas
derN = [dN1dl1 dN2dl1 dN3dl1 dN4dl1 dN5dl1 dN6dl1; ...
dN1dl2 dN2dl2 dN3dl2 dN4dl2 dN5dl2 dN6dl2];

% Calculo de la matriz de derivadas de las funciones de forma respecto a

% los ejes coordenados
dN = (inv(J1)*derN)';

% Calculo de la matriz B
B = [dN(1,1) 0 dN(2,1) 0 dN(3,1) 0 dN(4,1) 0 dN(5,1) 0 dN(6,1) 0; ...
0 dN(1,2) 0 dN(2,2) 0 dN(3,2) 0 dN(4,2) 0 dN(5,2) 0 dN(6,2); ...
dN(1,2) dN(1,1) dN(2,2) dN(2,1) dN(3,2) dN(3,1) dN(4,2) dN(4,1) dN(5,2)
dN(5,1) dN(6,2) dN(6,1)];
end
%-----------------------------------------------------------------------------

2.2 Integración numérica

La ventaja de utilizar coordenadas de área está en la capacidad de evaluar las

integrales, como en el caso de una dimensión. Para el elemento bidimensional con dos
coordenadas ( L1 , L2 ) se tiene

L
a!b!
∫ L L dx = (1 + a + b + c)! 2 A
a b
1 2 (2.16)
0

La integral de longitud L se define entre dos nodos a , b . esta relación también será
válida en la definición de las integrales que son sólo una función de la longitud a lo
largo de un borde de un elemento. Para los dos elementos tridimensionales con L1 ,
L2 , y L3 se tiene:

16
 L
a!b!
∫ L L L dA = (2 + a + b + c)! 2 A
a b c
1 2 3 (2.17)
0

Cuando A denota el área. Por ejemplo, en un triángulo lineal, el área de integración

de 1  2 esta dada por:

L L
1!1!0! A
∫ 1  2 dA = ∫ L1L2 L3dA 2A =
1 1 0
(2.18)
0 0
(2 + 1 + 1 + 0)! 12

Las integrales de área que se plantean son las siguientes

1 1− L

∫ ∫ f ( L , L ,L ) J dL L
0 0
1 2 3 1 2 (2.19)

Cuando las funciones de orden superior de interpolación (cuadrática, cúbica, entre

otras.) están escritas en función de las coordenadas de área, las matriz de las
integrales debe ser evaluada numéricamente, esto se debe a la inversa de la matriz
Jacobiana, que es una función racional propia de las coordenadas de área, y aparece
en la ecuación (1.22), cuando la derivada de las funciones de forma están
involucradas. Las formulas utilizadas en la ecuación (1.26) son válidas para las
integrales simples. Sin embargo, para integrales complejas es mejor dejar que el
computador haga la integración a fin de evitar los errores.

En el método de elementos finitos, se supone como:

1 1− L2
1 M
∫∫ f ( L1 ,L2 , L3 ) J dL1dL2 = ∑ wm g (( L1 ) m , ( L2 ) m , ( L3 ) m ,)
2 m =1
( 2.20 )
0 0

Donde g (⋅) Incluye J , el factor ½ corresponde al área en el sistema de coordenadas

local, ( L1 ) m denota los puntos específicos en el triangulo, y Wm son los pesos
asociados. El orden de la fórmula de cuadratura debe ser del menor número entero
mayor que la suma de las influencias de las coordenadas L1 , L2 , y L3 . Por ejemplo,
para integrar el producto L1 L2 L3 , la suma de los exponentes es tres, y se utiliza un
esquema de integración de cuarto grado.

17
Tabla1: Muestra los pesos y los valores de las coordenadas para las fórmulas de
integración numérica de diversos órdenes.

18
A continuación se presenta el código en Matlab par la construcción de la integración
numérica de un triangulo de 6 nodos.
%-----------------------------------------------------------------------------
%%Calcular k de del elemento
material_element = elements(element,end); % Numero de material
thickness_element = 1;
%Valores de L para la integración numerca
L =[1/2 0 1/2; 1/2 1/2 0; 0 1/2 1/2];
w = 1/3; %peso de Gauss
sum = zeros(12);
B_aux = zeros(3,12);
for p=1:3
[B detJ1] = calcular_B(x,y,L,p);
f = detJ1*B'*D{material_element}*B ;
sum = sum + f;
B_aux = B_aux + B;
end
k_element = 0.5*w*sum*thickness_element;

%-----------------------------------------------------------------------------

3. Elemento rectangular lineal bidimensional

3.1 Funciones de forma

En su forma más simple, el cuadrilátero se convierte en un elemento rectangular con

un sistema de coordenadas (en el espacio físico). Utilizando un sistema de
coordenadas locales, (naturales), tenemos un cuadrilátero generalizado.

La función de interpolación bilineal de un rectángulo de cuatro nodos esta dada como

φ = α 1 + α 2 x + α 3 y + α 4 xy (3.1)

El término xy fue elegido en lugar de x 2 o y 2 dado que xy no cambia el sistema y

supone una variación lineal de φ a lo largo de x o y . Para x constante, el elemento
es lineal en y , y viceversa, por lo tanto el nombre del elemento es bilineal.

∂φ
= α2 +α4 y
∂x
∂φ
= α3 + α4 x (3.2)
∂y

Fig.7 muestra el elemento de configuración rectangular y nodal. Si las distancias al

origen de coordenadas, que se encuentra en el punto medio del elemento, se denotan
por a y b , los valores de nodo se dan como

19
 .

Fig.7: Rectángulo de cuatro nodos

φ = φ1 para x = −b y = −a
φ = φ2 para x=b y = −a
(3.3)
φ = φ3 para x=b y=a
φ = φ4 para x = −b y=a

Sustituyendo estos valores en la ecuación (3.1) tenemos

φ1 + φ 2 + φ3 + φ 4
α1 =
4
φ1 − φ 2 − φ3 + φ 4
α2 =
4a
φ1 + φ 2 − φ3 − φ 4
α3 =
4b
φ1 − φ 2 + φ3 − φ 4
α4 = (3.4)
4ab

Si φ se aproxima en términos de los valores nodales y las funciones de forma como

se menciona anteriormente, se obtienen

φ = 1φ1 +  2φ 2 +  3φ3 +  4φ 4 (3.5)

Las funciones de forma son las siguientes

1
1 = (b − x)(a − y )
4ab
1
2 = (b + x)(a − y )
4ab
1
3 = (b + x)(a + y )
4ab
1
4 = (b − x)(a + y ) (3.6)
4ab

Para poder modelar estos elementos sin importar su ubicación y características se

procede a hacer una transformación que se presenta a continuación

20
3.2 Sistema de coordenadas naturales

En la mayoría de las situaciones, el dominio físico a ser modelado no consiste en los

costados rectos, fronteras ortogonales. Una región limitada por curvas puede ser
desacreditada con elementos cuadriláteros con lados curvos para obtener una
solución más precisa. La transformación de la recta de lados curvos se logra expresar
las coordenadas x , y , en términos de coordenadas curvilíneas:

x = x(r , s)
y = y (r , s) (3.6)

La elección de r y s depende de la geometría del elemento. El sistema de

coordenadas r , s se conoce normalmente como sistema de coordenadas naturales,
cuando el sistema de coordenadas se encuentran dentro del rango − 1 ≤ r ≤ 1 ,
− 1 ≤ s ≤ 1 . Para el caso rectangular, podemos encontrar el conjunto de dimensiones
coordenadas r , s tales que:

x
r=
b
y
s= (3.7)
a

Donde r y s efectivamente sustituirán la x y y que son las coordenadas globales

por coordenadas locales para un elemento individual, se había referido con
anterioridad a estos coeficientes específicos como índices de longitud (Fig. 7) por lo
tanto, la utilización de nuestra definición de la ecuación (3.7), se determinan las
funciones de forma bilineal como:

1
1 = (1 − r )(1 − s )
4
1
 2 = (1 + r )(1 − s )
4
1
 3 = (1 + r )(1 + s )
4
1
 4 = (1 − r )(1 + s ) (3.8)
4

El sistema de coordenadas local siempre varía entre -1 y 1. Tomamos ventaja de este

hecho cuando se formulen los términos derivados. Fig.8 muestra el sistema de
coordenadas r , s esta centrado en el elemento. Hay que tener en cuenta que no es
importante para este nuevo sistema de coordenadas que sean ortogonales también
puede ser curvilínea. Sin embargo, las derivadas con respecto a x y y no se pueden
obtener directamente y requieren de una transformación adicional. Vamos a tratar este
punto más adelante, por ahora, vamos a suponer un sistema simple de coordenadas
cartesianas.

21
 Fig.8: Sistema de coordenadas naturales

La evaluación de las derivadas ∂ i ∂x , ∂ i ∂y es similar a las derivadas de los

elementos triangulares mencionada anteriormente. La matriz Jacobiana puede ser
usada para encontrar las derivadas a través del inverso del Jacobiano con respecto a
la x , y y determinar la transformación.

Ejemplo 2: los valores de ∂ 1 ∂x y ∂ 1 ∂y en x = 2 y y = 2 para el elemento que se

muestra en la Fig.8. La transformación de x y y se expresan de la siguiente forma

1
x= (3s + 5)
2
y =r+2 (3.9)

Al evaluar la matriz Jacobiana obtenemos:

 ∂x ∂y 
 ∂s  = 1 3 0
J =  ∂s
∂y  2 0 2
(3.10)
∂x
 
 ∂r ∂r 

La inversa del Jacobiano esta dada por:

1 2 0
J −1 =   (3.11)
3  0 3

Las derivadas ∂ 1 ∂x y ∂ 1 ∂y se puede calcular una vez que la derivada ∂ 1 ∂s y

∂ 1 ∂r sean conocidas. Sabemos que

(1 − s)(1 − r )
1 = (3.12)
4

Las derivadas están dadas por

∂ 1 (1 − r )
=
∂s 4
∂ 1 (1 − s )
= (3.13)
∂r 4

22
Despejando de la ecuación (3.9) y remplazando x = 2 y y = 2 tenemos

1 1
s = ( 2 x − 5) = −
3 3
r = y−2=0 (3.14)

Así

∂ 1 (1 − 0) 1
= =−
∂s 4 4
 1
1 − 
∂ 1  3  1
= =− (3.15)
∂r 4 3

Remplazando obtenemos

 ∂ 1   1  1
 ∂x  1 2 0 −  − 
 ∂  =  4 6
0 3  ⋅  1 =  1 (3.16)
 1
 3   −  − 
 ∂y   3  3

A continuación se presenta el código en Matlab para la construcción de la matriz B

de un elemento cuadrilátero de 4 nodos
%-----------------------------------------------------------------------------
function [B detJ1] = calcular_B(r,s,x,y)

% Calculo de las derivadas de x e y respecto a los ejes coordenados

dxdr = 1/4*(1+s)*x(1)-1/4*(1+s)*x(2)-1/4*(1-s)*x(3)+1/4*(1-s)*x(4);
dxds = 1/4*(1+r)*x(1)+1/4*(1-r)*x(2)-1/4*(1-r)*x(3)-1/4*(1+r)*x(4);
dydr = 1/4*(1+s)*y(1)-1/4*(1+s)*y(2)-1/4*(1-s)*y(3)+1/4*(1-s)*y(4);
dyds = 1/4*(1+r)*y(1)+1/4*(1-r)*y(2)-1/4*(1-r)*y(3)-1/4*(1+r)*y(4);

% Calculo de las derivadas de las funciones de forma respecto a las

% coordenadas parametrizadas
dN1dr = 1/4 + 1/4*s;
dN1ds = 1/4 + 1/4*r;
dN2dr = -1/4 - 1/4*s;
dN2ds = 1/4 - 1/4*r;
dN3dr = -1/4 + 1/4*s;
dN3ds = -1/4 + 1/4*r;
dN4dr = 1/4 - 1/4*s;
dN4ds = -1/4 - 1/4*r;

% Almacenamiento del jacobiano

J1 = [dxdr dydr; ...
dxds dyds];

detJ1 = det(J1);

% Almacenamiento de la matriz de derivadas respecto a coordenadas

parametrizadas
derN = [dN1dr dN2dr dN3dr dN4dr; ...
dN1ds dN2ds dN3ds dN4ds];

23
 % Calculo de la matriz de derivadas de las funciones de forma respecto a
% los ejes coordenados
dN = (inv(J1)*derN)';

% Calculo de la matriz B
B = [dN(1,1) 0 dN(2,1) 0 dN(3,1) 0 dN(4,1) 0; ...
0 dN(1,2) 0 dN(2,2) 0 dN(3,2) 0 dN(4,2); ...
dN(1,2) dN(1,1) dN(2,2) dN(2,1) dN(3,2) dN(3,1) dN(4,2) dN(4,1)];

end
%-----------------------------------------------------------------------------

4. Elemento rectangular cuadrático bidimensional

4.1 Funciones de forma

Se puede desarrollar un modelo dos dimencional, de orden superior, para un elemento

cuadrilátero, tal como lo hicimos para el elemento triangular. Hay dos formas en que
esto se puede hacer que conducir a dos diferentes elementos de segundo grado: uno
con ocho nodos, uno en cada esquina y otra en el punto medio de cada lado, y el otro
cuenta con nueve nodos, siendo el noveno en el centro de la elemento este último es
sobre el que vamos a realizar la construcción ya que este se utilizan
predominantemente en la industria y se ha comprobado que es suficiente para la
mayoría de los efectos de cálculo.

El polinomio de interpolación para el elemento de ocho nodos, se muestra en la Fig.9

Consta de ocho términos.

La variable desconocida φ se expresa como

φ = α 1 + α 2 x + α 3 y + α 4 xy + α 5 x 2 + α 6 y 2 + α 7 x 2 y + α 8 xy 2 + ( 4.1)

Fig.9: Elemento cuadrilátero de ocho nodos

Si los valores del nodo se definen como en la Fig. 9.

24
 φ = φ1 para x = −b y = −a
φ = φ2 para x=0 y = −a
φ = φ3 para x=b y = −a
φ = φ4 para x=b y=0
φ = φ5 para x=b y=a
φ = φ6 para x=0 y=a
(4.2)
φ = φ7 para x = −b y=a
φ = φ8 para x = −b y=0

Sustituyendo los valores de la ecuación (2.8) a (2.7) los valores de las siguientes ocho
ecuaciones (expresado en forma de matriz):

 1 −b −a ab b2 a2 − ab 2 − a 2b α 1  φ1 
  α  φ 
 1 0 −a 0 0 a2 0 0   2  2
 1 b −a − ab b2 a2 − ab 2 − a 2b α 3  φ3 
     
 1 b 0 0 b2 0 0 0  α 4  = φ 4  (4.3)
 1 b a ab b2 a2 − ab 2 − a 2b α 5  φ5 
     
 1 0 a 0 0 a2 0 0  α 6  φ 6 
 1 −b a − ab b 2
a 2
− ab 2 − a 2b α  φ 
   7  7
 1 −b 0 0 b2 0 0 0  α 8  φ8 

Desde que el α i ' s con i = 1,...,8 , se puede obtener.

La aproximación de φ en función de los valores nodales y de cinco funciones de forma

φ =  1φ1 +  2φ 2 +  3φ3 +  4φ 4 +  5φ5 +  6φ 6 +  7φ 7 +  8φ8 ( 4 .4 )

Las funciones de forma son las siguientes

1
1 = (b − x)(a − y )(ab + ca + yb)
4(ab) 2
1
2 = (b 2 + x 2 )(a − y )
4ab 2
1
3 = (b + x)(a − y )(ax + by − ab)
4(ab) 2
1
4 = 2 (a 2 − y 2 )(b + x)
4a b
1
5 = (b + x)(a + y )(ax + by − ab)
4(ab) 2
1
6 = (b 2 − x 2 )(a + y )
4ab 2
1
7 = (b − x)(a + y )(ab + ax − by )
4(ab) 2
1
8 = 2 (a 2 − y 2 )(b − x) (4.5)
4a b

25
A continuación se presenta el código en Matlab para la construcción de la matriz B
de un elemento cuadrilátero de 8 nodos
%-----------------------------------------------------------------------------
function [B detJ1] = calcular_B(r,s,x,y)

% Calculo de las derivadas de x e y respecto a los ejes

% coordenados
dxdr = 1/4*((1+s)*x(1)-(1+s)*x(2)-(1-s)*x(3)+(1-s)*x(4)-2*(1+s)*(2*r)*x(5)
-2*(1-s^2)*x(6)-2*(1-s)*(2*r)*x(7)+2*(1-s^2)*x(8));
dxds = 1/4*((1+r)*x(1)+(1-r)*x(2)-(1-r)*x(3)-(1+r)*x(4)+2*(1-r^2)*x(5)-2*
(1-r)*(2*s)*x(6)-2*(1-r^2)*x(7)-2*(1+r)*(2*s)*x(8));
dydr = 1/4*((1+s)*y(1)-(1+s)*y(2)-(1-s)*y(3)+(1-s)*y(4)-2*(1+s)*(2*r)*x(5)
-2*(1-s^2)*x(6)-2*(1-s)*(2*r)*x(7)+2*(1-s^2)*x(8));
dyds = 1/4*((1+r)*y(1)+(1-r)*y(2)-(1-r)*y(3)-(1+r)*y(4)+2*(1-r^2)*x(5)-2*
(1-r)*(2*s)*x(6)-2*(1-r^2)*x(7)-2*(1+r)*(2*s)*x(8));

% Calculo de las derivadas de las funciones de forma respecto a las

% coordenadas parametrizadas
dN1dr = 1/4*(1+s);
dN1ds = 1/4*(1+r);
dN2dr = -1/4*(1+s);
dN2ds = 1/4*(1-r);
dN3dr = -1/4*(1-s);
dN3ds = -1/4*(1-r);
dN4dr = 1/4*(1-s);
dN4ds = -1/4*(1+r);
dN5dr = -1/2*(1+s)*(2*r);
dN5ds = 1/2*(1-r^2);
dN6dr = -1/2*(1-s^2);
dN6ds = -1/2*(1-r)*(2*s);
dN7dr = -1/2*(1-s)*(2*r);
dN7ds = -1/2*(1-r^2);
dN8dr = 1/2*(1-s^2);
dN8ds = -1/2*(1+r)*(2*s);

% Almacenamiento del jacobiano

J1 = [dxdr dydr; ...
dxds dyds];

detJ1 = det(J1);

% Almacenamiento de la matriz de derivadas respecto a coordenadas

parametrizadas
derN = [dN1dr dN2dr dN3dr dN4dr dN5dr dN6dr dN7dr dN8dr ; ...
dN1ds dN2ds dN3ds dN4ds dN5ds dN6ds dN7ds dN8ds ];

% Calculo de la matriz de derivadas de las funciones de forma respecto a

% los ejes coordenados
dN = (inv(J1)*derN)';

% Calculo de la matriz B
B=[dN(1,1) 0 dN(2,1) 0 dN(3,1 0 dN(4,1) 0 dN(5,1) 0 dN(6,1) 0 dN(7,1) 0
dN(8,1) 0; 0 dN(1,2) 0 dN(2,2) 0 dN(3,2) 0 dN(4,2) 0 dN(5,2) 0 dN(6,2)
0 dN(7,2) 0 dN(8,2);dN(1,2) dN(1,1) dN(2,2) dN(2,1) dN(3,2) dN(3,1)
dN(4,2) dN(4,1) dN(5,2) dN(5,1) dN(6,2) dN(6,1) dN(7,2) dN(7,1) dN(8,2)
dN(8,1)];

%-----------------------------------------------------------------------------

26
4.2 Integración numérica por medio de la cuadratura de Gauss

La formulación de las integrales dependen de r y s los límites de integración simple.

Esto es especialmente ventajoso cuando se trata de los límites de curvas o límites que
no son paralelos a los ejes de coordenadas. Las integrales se definen como
a b 1 1

∫ ∫ F ( x, y)dxdy = ∫ ∫ f (r, s) J drds

− a −b −1 −1
(4.6)

La matriz se convierte así en

1 1
 ∂ i ∂ j ∂ i ∂ j 
K= ∫ ∫ K  ∂x
−1 −1
∂x
+
∂y ∂y 
 dxdy

1 1
 ∂ i ∂ j ∂ i ∂ j 
= ∫ ∫ K  ∂x
−1 −1
∂x
+
∂y ∂y 
 J drds ( 4 .7 )

Donde las derivadas ∂ i ∂x y ∂ i ∂y se pueden obtener de la ecuación 3.10 en

función de r , s , ∂ i ∂r y ∂ i ∂s . Todas las integrales de área se convierten el
producto  i  j como se muestra a

1 1

∫∫ Ω
 i  j dxdy = ∫∫
−1 −1
i  j J drds ( 4 .8)

Integración del termino derecho esta dado por la ecuación 4.7, ya no implica
polinomios simples debido a la condición 1 / J de la expresión de las derivadas con
respecto a r y s , puede llegar a ser molesto. La práctica más común es utilizar el
procedimiento de integración numérica se conoce como la cuadratura de Gauss. El
procedimiento es fácil de programar y tiene soluciones muy aproximadas para las
integrales

Una cuadratura de Gauss aproxima una integral definida por medio de una suma
ponderada de un conjunto finito de puntos. Por ejemplo, en una dimensión tenemos
1 

∫ f ( r ) dr ≅ ∑W
i =1
i f ( ri ) ( 4 .9 )
−1

Donde  es el número de puntos de integración, Wi son los pesos de Gauss, y los

Ri son los puntos de Gauss. Con los puntos de Gauss, podemos integrar
exactamente un polinomio de grado 2  − 1 . Así pues, si f (r ) es un polinomio
cúbico f (r ) = a + br + cr 2 + dr 3 , entonces la elección  = 2 , los valores de Ri se
definen como

27
 r1 = −1 / 3
r2 = 1 / 3 ( 4.10 )

Las ponderaciones serian W1 = W2 = 1 . Para las integrales dobles en la ecuación 4.8,

tenemos.
1 1  

∫∫ f ( r , s ) drds ≅ ∑ ∑ WiW j f ( ri , s j ) ( 4.11)

−1 −1 i =1 j =1

Los puntos de Gauss están definidos por los siguientes valores

r1 = −1 / 3
r2 = 1 / 3
s1 = −1 / 3
s2 = 1 / 3 ( 4.12 )

Para cuatro puntos de Gauss se utilizan los cuatro elementos del cuadrilátero como se
muestra en la Fig. 10 y el tabla2

Fig.10: Cuadrilátero con cuatro nodos y cuatro puntos de Gauss

Tabla2: Integración con cuatro puntos de Gauss

Si el elemento cuadrático es de ocho nodos o de nueve se utilizan la integral que

aparece en la ecuación 4.7 se quedando definida como :
1 1 3 3

∫ ∫ f (r , s )drds ≅ ∑ ∑ W W
−1 −1 i =1 j =1
i j f ( ri , s j ) ( 4.13)

28
Los nueve puntos de integración son presentados en la Fig11 y los pesos se muestran
en la tabla 3

Fig.11: Integración numérica con nueve puntos de Gauss

Tabla3: Pesos de Gauss.

La solución obtenida con cuadriláteros bilineales implica una matriz 4x4, mientras que
un cuadrilátero cuadrático implica una matriz de 8x8. Además, los elementos bilineales
requieren una evaluación en 4 puntos de Gauss, y no involucran ecuaciones
cuadráticas, por lo que es nueve veces menos costoso en términos de cálculo para
construir la matriz de rigidez de elemento por elemento. Esto puede ser una
consideración importante al escoger el tipo de elemento que se utilizará, ya que
significa una diferencia de casi un orden de magnitud en el tiempo computacional.

A continuación se presenta el código en Matlab para la construcción de la integración

numérica para un cuadrilátero de 4 nodos
%-----------------------------------------------------------------------------
pg = [sqrt(1/3); -sqrt(1/3)]; %pesos de gauss

%%Calcular k de del elemento

material_element = elements(element,end);
thickness_element = 1;
%matris k de un elemento
aux = 0;
for ri = 1:2
r = pg(ri);

29
 for si = 1:2
s = pg(si);
%%Calcular B del elemento
[B detJ1] = calcular_B(r,s,x,y);
k_element = k_element +
B'*D{material_element}*B*detJ1*thickness_element;
aux = aux + B;
end
end
%-----------------------------------------------------------------------------

A continuación se presenta el código en Matlab para la construcción de la integración

numérica para un cuadrilátero de 8 nodos
%-----------------------------------------------------------------------------
pg = [0.77459669241483; 0 ; -0.77459669241483];

%%Calcular k de del elemento

material_element = elements(element,end);
thickness_element = 1;
k_element = zeros(16,16);
%matris k de un elemento
aux = 0;
for ri = 1:3
r = pg(ri);
for si = 1:3
s = pg(si);
%%Calcular B del elemento
[B detJ1] = calcular_B(r,s,x,y);
k_element = k_element B'*D{material_element}*B*detJ1*thickness_element;
aux = aux + B;
end
end
%-----------------------------------------------------------------------------

5. Elementos tridimensionales

Los elementos más comunes tridimensionales provienen de variaciones de los

elementos bidimensionales, presentados anteriormente, es decir, el tetraedro es un
elemento lineal limitado por los lados rectos; pero el hexaedro que es de orden
superior puede tener superficies curvas. Esto mismo ocurre con los elementos
tetraédricos y los cilindros.

Fig.12: Algunos elementos tridimensionales

30
6. Tetraedro lineal

Un tetraedro lineal esta definido por cuatro nodos como se muestra en la Fig.13, y su
función de interpolación esta dada por:

φ = α1 + α 2 x + α 3 y + α 4 z ( 6 . 1)

Fig.13: Elemento Tetraedro

Para los cuatro nodos situados en los vértices del tetraedro, se pude escribir el
siguiente sistema de ecuaciones

φ1 = α 1 + α 2 x1 + α 3 y 1 + α 4 z 1
φ2 = α 1 + α 2 x2 + α 3 y2 + α 4 z2
φ3 = α 1 + α 2 x3 + α 3 y3 + α 4 z3
φ4 = α 1 + α 2 x4 + α 3 y4 + α 4 z4 ( 6 .2 )

de forma matricial tenemos

φ = Cα ( 6 .3)

Que también se puede escribir como

 φ 1  1 x1 y1 z1   α 1 
φ  1 x2 y2 z 2  α 2 
 2 =  ⋅ ( 6 .4 )
φ 3  1 x3 y3 z 3  α 3 
     
φ 4  1 x4 y4 z 4  α 4 

Tomando la inversa de C , el vector de la columna α se obtiene como

α 1 
α 
φ = [1 x y z ] ⋅  2  = [1 x y z ] ⋅ C −1 φ ( 6 .5 )
α 3 
 
α 4 

31
La cantidad [1 x y z ] ⋅ C − 1 es en realidad la forma de las funciones matriciales, es
decir,

φ = φ ( 6 .6 )

Donde

 = [1 x y z ]⋅ C −1
( 6 .7 )

Un sistema de coordenadas de volumen para el tetraedro se puede establecer de la

misma manera como el sistema coordenado de áreas para un triangulo de dos
dimensiones, en el hexaedro de tres dimensiones, requiere cuatro relaciones de
distancia, cada una normal a un lado: L1 , L 2 , L 3 y L 4 . El volumen V1 se muestra en
la Fig.14.

Fig.14: Sistema coordenado de volumen

Las cuatro coordenadas están relacionadas de la siguiente forma

L1 + L 2 + L 3 + L 4 = 1 ( 6 .8 )

La forma lineal de la función queda definida por

 1 = L1
 2 = L2
 3 = L3
 4 = L4 ( 6 .9 )

Las integrales de volumen pueden ser evaluadas fácilmente de la siguiente relación

a ! b! c! d !
∫L Lb2 Lc3 L d4 dV =
a
1 6V ( 6 . 10 )
V
( 3 + a + b + c + d )!

Del mismo modo que se utilizo para los elementos triangulares.

32
 Fig.15: Tetraedro cuadrático

La funciones de forma de los elementos de orden superior son similares a las del
tetraedro lineal, Por ejemplo, la función de forma en el nodo 1 para el tetraedro de
segundo grado, que se muestra en la Fig.15 queda definida por

 1 = ( 2 L 1 − 1) L 1 ( 6 . 11 )

Cuando L1 = 0 pasa a través de los nodos 2, 3, 4, 6, 9 y 10. Para L1 = 1 / 2 , el plano

pasa a través de los nodos 5, 7 y 8. En el nodo 5,

 5 = 4 L1 L 2 ( 6 . 12 )

Cuando L 2 = 0 pasa a través de los nodos 1, 3, 4, 7, 8 y 10. La función de las

relaciones de forma adicional puede ser fácilmente evaluada por referencia a la
ecuación 4.26. Que da las funciones de forma para el elemento triangular de dos
dimensiones de segundo grado.

Ejemplo 3: Determinar las funciones de forma para el elemento de tetraedro lineal que
se muestra en la Fig. 16 de la ecuación 6.4.

Fig.16: Tetraedro lineal

33
En primer lugar evaluaremos C :

1 x1 y1 z 1  1 3 2 0
1 x2 y2 z 2  1 0 4 0 
C =  = ( 6 . 13 )
1 x3 y3 z 3  1 0 0 0
   
1 x4 y4 z 4  1 1 2 5

Se hace necesario encontrar la inversa C − 1 a fin de obtener las funciones de forma

como se define por la ecuación 6.7. La inversa se obtiene de la siguiente manera

 0 0 60 0 
 20 − 10 − 10 0 
 
 0 15 − 15 0 
 
− 4 − 4 − 4 12
=  
adjC
C = ( 6 . 14 )
c 60

Las funciones de forma son dadas por el siguiente producto

 0 0 60 0
 0 
 = [1 x y z ] ⋅ C −1 =
[1 x y z ]  20 − 10 − 10
60  0 15 − 15 0
 
− 4 −4 −4 12 

1
= [( 20 x − 4 z )( − 10 x + 15 y − 4 z )( 60 − 10 x − 15 y − 4 z )(12 z ) ] ( 6 .15 )
60

Así, quedan definidas las funciones de forma para cada nodo

1
1 = (5 x − z )
15
1
2 = ( − 10 x + 15 y − 4 z )
60
1
3 = ( 60 − 10 x − 15 y − 4 z )
60
z
4 = ( 6 . 16 )
5

A continuación se presenta el código en Matlab para la construcción la matriz B y el

volumen V para un tetra de 4 nodos

%-----------------------------------------------------------------------------
%calcular_V
function V = calcular_V(coordenate_node1,coordenate_node2,coordenate_node3,
coordenate_node4)

m = [coordenate_node1-coordenate_node2,coordenate_node2- coordenate_node3,
coordenate_node3-coordenate_node4];

V = (1/6)*abs(det(m));

34
end

%-----------------------------------------------------------------------------
function B = calcular_B(coordenate_node1,coordenate_node2,coordenate_node3,
coordenate_node4)

% Calculo de las derivadas de de x e y respecto a los ejes coordenados

x(1) = coordenate_node1(1);
x(2) = coordenate_node2(1);
x(3) = coordenate_node3(1);
x(4) = coordenate_node4(1);
y(1) = coordenate_node1(2);
y(2) = coordenate_node2(2);
y(3) = coordenate_node3(2);
y(4) = coordenate_node4(2);
z(1) = coordenate_node1(3);
z(2) = coordenate_node2(3);
z(3) = coordenate_node3(3);
z(4) = coordenate_node4(3);
dxdr = -1*x(1)+x(2);
dxds = -1*x(1)+x(3);
dxdt = -1*x(1)+x(4);
dydr = -1*y(1)+y(2);
dyds = -1*y(1)+y(3);
dydt = -1*y(1)+y(4);
dzdr = -1*z(1)+z(2);
dzds = -1*z(1)+z(3);
dzdt = -1*z(1)+z(4);

% Calculo de las derivadas de las funciones de forma respecto a las

% coordenadas parametrizadas
dN1dr =-1;
dN1ds =-1;
dN1dt =-1;
dN2dr = 1;
dN2ds = 0;
dN2dt = 0;
dN3dr = 0;
dN3ds = 1;
dN3dt = 0;
dN4dr = 0;
dN4ds = 0;
dN4dt = 1;

% Almacenamiento del jacobiano

J1 = [dxdr dydr dzdr; dxds dyds dzds; dxdt dydt dzdt];
% Almacenamiento de la matriz de derivadas respecto a coordenadas
parametrizadas
derN = [dN1dr dN2dr dN3dr dN4dr;
dN1ds dN2ds dN3ds dN4ds;
dN1dt dN2dt dN3dt dN4dt];

% Calculo de la matriz de derivadas de las funciones de forma respecto a

% los ejes coordenados
dN = inv(J1)*derN;

B = [dN(1,1) 0 0 dN(1,2) 0 0 dN(1,3) 0 0 dN(1,4) 0 0; 0 dN(2,1) 0 0 dN(2,2)

0 0 dN(2,3) 0 0 dN(2,4) 0; 0 0 dN(3,1) 0 0 dN(3,2) 0 0 dN(3,3) 0 0
dN(3,4); dN(2,1) dN(1,1) 0 dN(2,2) dN(1,2) 0 dN(2,3) dN(1,3) 0 dN(2,4)
dN(1,4) 0; 0 dN(3,1) dN(2,1) 0 dN(3,2) dN(2,2) 0 dN(3,3) dN(2,3) 0
dN(3,4) dN(2,4); dN(3,1) 0 dN(1,1) dN(3,2) 0 dN(1,2) dN(3,3) 0 dN(1,3)
dN(3,4) 0 dN(1,4)];
end
%-----------------------------------------------------------------------------

35
7. Hexaedro lineal

El tetraedro es un elemento simple, que se extiende directamente desde el elemento

triangular. Sin embargo, el tetraedro es difícil de visualizar en una compleja maraña
tridimensional y puede requerir un esfuerzo considerable para definir el elemento de
conectividad. Una solución más precisa con menos esfuerzo, pero más puntos de
integración, se puede lograr con el elemento de hexaedro. El elemento más simple de
esta familia es un hexaedro de ocho nodos, que se muestra en la Fig. 17

Fig.17: Hexaedro de ocho nodos

La función de interpolación esta dada por

φ = α 1 + α 2 x + α 3 y + α 4 z + α 5 xy + α 6 xz + α 7 yz + α 8 xyz ( 7 . 1)

La ecuación 7.1 se puede expresar en forma de matriz como se presenta a

continuación

 φ 1  1 x1 y1 z1 x1 y 1 x1 z 1 y1 z1 x1 y 1 z 1   α 1 
φ  1 x2 y2 z2 x2 y2 x2 z2 y2 z2 x 2 y 2 z 2  α 2 
 2 =  ⋅ ( 7 .2 )

 







 

     
 φ 8  1 x8 y8 z8 x8 y 8 x8 z 8 y8 z8 x 8 y 8 z 8  α 8 

Dado que α = C − 1φ , podemos encontrar los cuatro valores de la expresión

φ = [1 x y z xy xz yz xyz ] ⋅ C − 1 φ ( 7 .3)

A partir de qué forma las funciones se obtienen como se mostró anteriormente

 = [1 x y z xy xz yz xyz ] ⋅ C − 1 ( 7 .3)

Como en los elementos cuadriláteros bidimensionales, un sistema de coordenadas

naturales también puede ser diseñado para el hexaedro, como se muestra en la Fig.
18. Las funciones de forma quedan definidas así.

36
 Fig.18: Sistema de coordenadas naturales

 1   (1 − r )(1 − s )(1 − t ) 
   (1 + r )(1 − s )(1 − t ) 
 2 
3  (1 + r )(1 + s )(1 − t ) 
   
  4  = 1  (1 − r )(1 + s )(1 − t ) 
( 7 .4 )
  5  8  (1 − r )(1 − s )(1 + t ) 
   
6   (1 + r )(1 − s )(1 + t ) 
   (1 + r )(1 + s )(1 + t ) 
 7  
  8   (1 − r )(1 + s )(1 + t ) 

La ecuación 7,4 se puede escribir en forma más compacta como

1
i = (1 − rri )(1 − ss i )(1 − tt i ) con i = 1, 2 , … ,8 ( 7 .5 )
8

Cuando ri , s i y t i = ± 1 , dependiendo de la ubicación nodal.

Ejemplo 4: Determinar la expresión  1 con coordenadas r , s , t . De la ecuación 7.5.

1
1 = (1 − rr i )(1 − ss i )(1 − tt i ) ( 7 .6 )
8

Podemos ver de la a la Fig.15 que

ri = − 1 si = − 1 ti = − 1

Así,

1
1 = (1 − r )(1 − s )(1 − t ) ( 7 .7 )
8

La extensión de los 8 nodos de elemento de segundo grado a los 20 nodos hexaedro

cuadrática puede hacerse fácilmente mediante la adición de un nodo en la mitad de
cada borde del elemento. Esto puede verse en la Fig.19

37
 Fig.19: Hexaedro cuadrático

Las funciones de forma para los nodos de esquina se obtienen de la siguiente relación

1
i = (1 − rri )(1 − ss i )(1 − tt i )( rri + ss i + tt i − 2 ) con i = 1, 2 , … ,8 ( 7 .8 )
8

Para los nodos intermedios se tienen estas relaciones,

1
i = (1 − r 2 )(1 − ss i )(1 − tt i ) con i = 10 ,12 ,18 , 20 ( 7 .9 )
4
1
 i = (1 − rr i )(1 − s 2 )(1 − tt i ) con i = 9 ,11 ,17 ,19 ( 7 . 10 )
8
1
 i = (1 − rr i )(1 − ss i )(1 − t 2 ) con i = 13 ,14 ,15 ,16 ( 7 . 11 )
8

Para un sistema de coordenadas locales orientado como se muestra en la Fig.15. El

elemento hexaedro de 9 nodos se extiende a uno de 27 nodos con un nodo en el
centro de cada una de las seis caras y uno en el centro del elemento.

Observe que cuando una de las coordenadas natural es ± 1 , los tres elementos de
dimensiones se reducen a sus dos dimensiones homólogas. Asimismo, cuando dos de
las coordenadas de la iguales ± 1 , la forma de la función se reduce a la función de una
forma tridimensional, permitiendo tres elementos de dimensiones que se unió a los
elementos de línea. La continuidad entre los elementos de existir más de una cara de
elemento (área) de un elemento tridimensional, de los valores nodales describir la
variación de la variable (φ ) desconocida de forma idéntica en la cara elemento común
a los elementos adyacentes.

A continuación se presenta el código en Matlab para la construcción la matriz B para

un hexaedro de 8 nodos
%-----------------------------------------------------------------------------
function [B detJ1] = calcular_B(r,s,t,x,y)

% Calculo de las derivadas de x e y respecto a los ejes coordenados

dxdr = (1/2)*(x(1)-x(4));
dxds = 0;
dxdt = 0;
dydr = 0;

38
 dyds = (1/2)*(y(2)-y(1));
dydt = 0;
dzdr = 0;
dzds = 0;
dzdt = (1/2)*(z(5)-z(1));

% Calculo de las derivadas de las funciones de forma respecto a las

% coordenadas parametrizadas
dN1dr = (-1/8)*(1-s)*(1-t);
dN1ds = (-1/8)*(1-r)*(1-t);
dN1dt = (-1/8)*(1-r)*(1-s);
dN2dr = (1/8)*(1-s)*(1-t);
dN2ds = (-1/8)*(1+r)*(1-t);
dN2dt = (-1/8)*(1+r)*(1-s);
dN3dr = (1/8)*(1+s)*(1-t);
dN3ds = (1/8)*(1+r)*(1-t);
dN3dt = (-1/8)*(1+r)*(1+s);
dN4dr = (-1/8)*(1+s)*(1-t);
dN4ds = (1/8)*(1-r)*(1-t);
dN4dt = (-1/8)*(1-r)*(1+s);
dN5dr = (-1/8)*(1-s)*(1+t);
dN5ds = (-1/8)*(1-r)*(1+t);
dN5dt = (1/8)*(1-r)*(1-s);
dN6dr = (1/8)*(1-s)*(1+t);
dN6ds = (-1/8)*(1+r)*(1+t);
dN6dt = (1/8)*(1+r)*(1-s);
dN7dr = (1/8)*(1+s)*(1+t);
dN7ds = (1/8)*(1+r)*(1+t);
dN7dt = (1/8)*(1+r)*(1+s);
dN8dr = (-1/8)*(1+s)*(1+t);
dN8ds = (1/8)*(1-r)*(1+t);
dN8dt = (1/8)*(1-r)*(1+s);

% Almacenamiento del jacobiano

J1 = [dxdr dydr dzdr; ...
dxds dyds dzds;...
dxdt dydt dzdt];

detJ1 = det(J1);

% Almacenamiento de la matriz de derivadas respecto a coordenadas

parametrizadas
derN = [dN1dr dN2dr dN3dr dN4dr dN5dr dN6dr dN7dr dN8dr; ...
dN1ds dN2ds dN3ds dN4ds dN5ds dN6ds dN7ds dN8ds;...
dN1dt dN2dt dN3dt dN4dt dN5dt dN6dt dN7dt dN8dt];

% Calculo de la matriz de derivadas de las funciones de forma respecto a

% los ejes coordenados
dN = (inv(J1)*derN)';

% Calculo de la matriz B
B = [dN(1,1) 0 0 dN(2,1) 0 0 dN(3,1) 0 0 dN(4,1) 0 0 dN(5,1) 0 0 dN(6,1) 0
0 dN(7,1) 0 0 dN(8,1) 0 0; 0 dN(1,2) 0 0 dN(2,2) 0 0 dN(3,2) 0 0
dN(4,2) 0 0 dN(5,2) 0 0 dN(6,2) 0 0 dN(7,2) 0 0 dN(8,2) 0; 0 0
dN(1,3) 0 0 dN(2,3) 0 0 dN(3,3) 0 0 dN(4,2) 0 0 dN(5,3) 0 0 dN(6,3) 0
0 dN(7,3) 0 0 dN(8,3); dN(1,2) dN(1,1) 0 dN(2,2) dN(2,1) 0 dN(3,2)
dN(3,1) 0 dN(4,2) dN(4,1) 0 dN(5,2) dN(5,1) 0 dN(6,2) dN(6,1) 0 dN(7,2)
dN(7,1) 0 dN(8,2) dN(8,1) 0; dN(1,3) 0 dN(1,1) dN(2,3) 0 dN(2,1) dN(3,3)
0 dN(3,1) dN(4,3) 0 dN(4,1) dN(5,3) 0 dN(5,1) dN(6,3) 0 dN(6,1) dN(7,3)
0 dN(7,1) dN(8,3) 0dN(8,1); 0dN(1,3) dN(1,2) 0dN(2,3) dN(2,2) 0 dN(3,3)
dN(3,2) 0 dN(4,3) dN(4,2) 0 dN(5,3) dN(5,2) 0 dN(6,3) dN(6,2) 0 dN(7,3)
dN(7,2) 0 dN(8,3) dN(8,2)];

end
%-----------------------------------------------------------------------------

39
8. Integración numérica

En este punto, la dificultad en la evaluación de las integrales de funciones de forma

relacionadas con más de unos pocos nodos que debería ser obvio. Utilizando las
coordenadas de volumen, los elementos tetraédricos se pueden evaluar exactamente
usando la ecuación. 6,10, sin embargo, la cantidad de álgebra asociados con la
multiplicación del producto interno se hace muy laboriosa y no rápida.
Las integrales de área que se producen en dos elementos triangulares son
dimensiones de la forma.

1 1− L 2
I = ∫∫ 0 0
f ( L 1 , L 2 , L 3 ) J dL 1 dL 2

Cuando las funciones de interpolación se expresan en términos de coordenadas de

área. Dado que el jacobiano esta en función de las coordenadas de área, una forma
explícita de la inversa del jacobiano es casi imposible. Si bien algunas integrales
pueden ser evaluadas usando la ecuación 2.17, que es las dos dimensiones de la
forma de la ecuación 7.14, la probabilidad de error es grande ya que el número de
términos aumenta dramáticamente.
n
1 1− L 2 1
∫∫ f ( L 1 , L 2 , L 3 ) J dL 1 dL 2 = ∑ w i g [( L 1 ) i , ( L 2 ) i , ( L 3 ) i ]
0 0 2 i =1

Si el orden de la integración ( n ) se determina mediante la suma de la potencia de las

coordenadas L1 , L2 y. L3

De manera similar, los elementos del tetraedro son evaluados como

1 1− L1 1− L 2
I = ∫∫
0 0 ∫
0
f ( L1 , L 2 , L 3 , L 4 ) J dL 1 dL 2 dL 3

n
1
=
6
∑ w g [( L ) , ( L
i =1
i 1 i 2 ) i , ( L3 ) i , ( L 4 ) i ] ( 8 . 1)

La integración numérica de los puntos nodales del tetraedros se muestran en la tabla4

es que queremos integrar el producto L1 L 2 L 4 , la suma de los exponenciales es 3, y
requiere una integración de cuarto orden ( n = 4 puntos ) .

40
 Tabla4: Puntos nodales para un tetraedro.

El procedimiento para la evaluación de las integrales del elemento para hexaedros es

casi idéntico al procedimiento utilizado para dos cuadriláteros tridimensionales. La
matriz jacobiana se convierte ahora en una matriz de 3x3, a diferencia de de la anterior
que era de 2x2 para el cuadrilátero. La forma general de la integración esta dada por
n n n

∑∑∑ww
1 1 1
∫ ∫ ∫
−1 −1 −1
f ( r , s , t ) J drds dt =
i =1 j =1 k =1
i j w k f ( ri , s j , t k ) J ( ri , s j , t k ) (8 . 2 )

w i , w j , wk son los factores de peso asociada con los puntos de Gauss respectivos.
Obsérvese que la ecuación. 8.2 se obtiene combinando una fórmula tridimensional
establecida originalmente para el único elemento tridimensional. Los factores de
ponderación y el punto de integración de las coordenadas se muestran en la tabla5.

Tabla5: Factores de ponderación y puntos de integración.

A continuación se presenta el código en Matlab para la construcción de la integración

numérica de un hexaedro de 8 nodos.

41
%-----------------------------------------------------------------------------
pg = [sqrt(1/3); -sqrt(1/3)]; %pesos de gauss

%%Calcular k de del elemento

material_element = elements(element,4);
thickness_element = 1;
k = zeros(24,24); %matris k de un elemento
for ri = 1:2
r = pg(ri);
for si = 1:2
s = pg(si);
for ti = 1:2
t = pg(ti);
%%Calcular B del elemento
[B detJ1] = calcular_B(r,s,t,x,y);
k = k + B'*D{material_element}*B*detJ1*thickness_element;
B_T{element} = B_T{element} + B; % Matrices B de cada elemento
end
end
end
%-----------------------------------------------------------------------------

9. Ensamble de la matriz de rigidez (K)

Para construir la matriz de rigidez (K ) es necesario construir la k de cada elemento ya

que al ensamblar cada una de estas matrices obtenemos la matriz de rigidez global.

La matriz k para cada elemento esta dada por:

k = ∫ B i′DB i d Ω ( 9 .1)
Ω

Donde B es una matriz que depende de cada tipo de elemento ya que esta definida
por las derivadas parciales de las funciones de forma con respecto a x y y (las
funciones de forma cambian para cada tipo de elemento como se presento
anteriormente), la matriz B queda definida por:

 ∂ i 
 0 
 ∂x 
∂ i 
Bi =  0 (9 .2 )
 ∂y 
 ∂ i ∂ i 
 
 ∂ y ∂ x 

Y D es la matriz de constantes elásticas definida a partir del modulo de Young (E ) y

el coeficiente de Poisson (v ) .

Matriz D para elementos planos:

1 v 0 
E  
D= v 1 0 ( 9 . 3)
1− v2  
 0 0 0 .5 (1 − v ) 

42
Matriz D para elementos de volumen:

 v v 
 1 1− v 1− v
0 0 0 
 v v 
 1 0 0 0 
1 − v 1− v 
 v v
1 0 0 0 
E (1 − v )  1 − v 1− v 
D=  1 − 2v  (9 .4 )
(1 + v )(1 − 2 v )  0 0 0 0 0 
 2 (1 − v ) 
 1 − 2v 
 0 0 0 0
2 (1 − v )
0 
 
 0 1 − 2v 
0 0 0 0
 2 (1 − v ) 

Para explicar el ensamblaje de cada k por elemento utilizaremos el siguiente ejemplo

Ejemplo 5: A partir de la Fig.20 construir la matriz de rigidez global K

Fig.20: Ejemplo 5 donde R indica una restricción y P un peso

Para el primer elemento que tiene los nodos 1, 2 construimos la matriz k

multiplicamos por el vector de desplazamientos y obtenemos el vector de restricciones
como se presenta en le ecuación (9.6).

EA
k1 = (9 .5 )
l1

 k1(1,1) k1(1, 2 )   u1   R 
k ⋅ =
k1( 2 , 2 )  u 2   0 
(9 .6 )
 1( 2 ,1)

Para el segundo elemento compuesto por los nodos 2 y 3 tenemos:

EA
k2 = (9 .7 )
l2

 k 2 (1,1) k 2 (1, 2 )  u 2   0 
k ⋅ =
k 2 ( 2 , 2 )   u 3   0 
(9 .8 )
 2 ( 2 ,1)

43
Para el tercer elemento compuesto por los nodos 3 y 4 tenemos:

EA
k3 = (9 .9 )
l3

 k 3 (1,1) k 3 (1, 2 )   u 3   0 
k ⋅ =
k 3 ( 2 , 2 )  u 4   P 
(9 .10 )
 3 ( 2 ,1)

Después de tener cada una de las k por cada elemento comenzamos a ensamblar la
matriz de rigidez global para esto construimos una matriz cuadrada del numero total de
desplazamientos, para este caso es una matriz de 4x4; la forma como se ensambla es
colocando la matriz k de tal forma que al multiplicar por el vector de desplazamientos
obtengamos el mismo sistema de ecuaciones, como se muestra a continuación

 k1(1,1) k1(1, 2 ) 0 0   u1   R 
k
 1( 2 ,1) k 1( 2 , 2 ) 0 0  u 2   0 
⋅ = (9 .11)
 0 0 0 0  u 3   0 
     
 0 0 0 0  u 4   P 

Esto se realiza para cada una de las matrices k , teniendo en cuenta que si hay algo
en la posición que le corresponde se suma a lo que había, como se muestra en la
ecuación (9.12)

 k1(1,1) k1(1, 2 ) 0 0   u1  R

k k1( 2 , 2 ) + k 2 (1,1) k 2 (1, 2 ) 0  u 2  0
 1( 2 ,1) ⋅  =   (9 .12 )
 0 k 2 ( 2 ,1) k 2 ( 2 , 2 ) + k 3 (1,1) k 3 (1, 2 )   u 3  0
     
 0 0 k 3 ( 2 ,1) k 3 ( 2 , 2 )  u 4  P

Otra forma de construir la matriz de rigidez global de manera más simple para trabajar
con muchos elementos será explicada a partir del ejemplo 6

Ejemplo 6: A partir de la Fig.21 construir la matriz de rigidez global K y encontrar la

solución para KU = F , donde U es el vector de deformaciones y F es el vector de
fuerzas.

Lo primero que vamos a realizar es una tabla que esta compuesta por el nodo, la
posición del nodo y las restricciones en x , y que este posee.

44
 Fig.21: Viga con 12 elementos triangulares

Tabla6: Presenta la posición y las restricciones que tiene cada nodo.

Después de esto se construye una tabla que nos permita saber que nodos pertenecen
a un mismo elemento, además de esto construimos una tabla que nos permita
identificar donde están las fuerzas a las que es sometida la pieza.

Tabla7: Nodos por elemento.

45
 Tabla8: Fuerzas por elemento

Para determinar la magnitud de la matriz de rigidez global es necesario construir una

tabla a partir de las restricciones que se mostraron en la Tabla6. Lo que se hace es
cambiar por los unos de las restricciones por ceros y los ceros de las restricciones
cambiarlos por el número de ceros que tenga en la tabla de las restricciones hasta esa
posición, como se muestra a continuación:

Tabla9: Matriz ID

El numero de ceros de la tabla de restricciones es 18 son los grados de libertad y por

esto la matriz de rigidez global es de 18x18.

Construimos ahora un vector LM que consiste en tomar la información perteneciente al

elemento de la matriz ID, si tomamos el primer elemento que tiene los nodos 1, 2 y 3
entonces tomamos los valores de la tabla9 en las filas 1,2 y 3, entonces el vector LM
queda definido como:

LM = [0 0 1 2 0 0] (9 .13 )

Después de esto se multiplica cada casilla por el vector LM y las coordenadas de los
productos que den diferentes de cero son las coordenadas del valor de la matriz k del
elemento que tiene que ser ensamblado en las posiciones del vector LM en la matriz
de rigidez global, para el caso del primer elemento al multiplicar la casilla tres del

46
vector LM por la casilla tres LM da diferente de cero, lo que indica que el valor en la
posición (3,3) de la matriz k del elemento es el valor que corresponde a la posición
(1,1) en la matriz de rigidez global.

Fig.22: Ensamble de de la matriz global K

Cuando dos valores tienen que ocupar la misma posición en la matriz de rigidez global
K , lo que se hace es sumar estos elementos.

A continuación se presenta el código en Matlab para el ensamble de la matriz de la

matriz de rigidez global K (para un elemento hexaedro)
NOTA: para los otros elementos solo es necesario cambiar el número de nodos.
%-----------------------------------------------------------------------------
% Calculo de ID y de grados de libertad libres
function [ID n_fdof] = calcular_ID(constraints)
[n_elements n_dof] = size(constraints);
ID = zeros(n_elements,n_dof);
n_fdof = 0;
for i = 1:n_elements
for j = 1:n_dof
if not(constraints(i,j))
n_fdof = n_fdof + 1;
ID(i,j) = n_fdof;
end
end
end
end

%%----------------------------------------------------------------------------

%%Ensamblaje de K
LM_element = [ ID(node1,:) ID(node2,:) ID(node3,:) ID(node4,:) ID(node5,:)
ID(node6,:) ID(node7,:) ID(node8,:)];
assemble_K(k_element,LM_element);

%%----------------------------------------------------------------------------

%%Nested functions:

%Ensamblaje de K
function assemble_K(k,LM)
%Guarda en m el número de columnas de la matriz k
m = length(k);
%Guarda el número de elementos del vector K_VAL (que es igual que K_COL y
K_ROW)
LENGTH_K = length(K_VAL);
for i=1:m

47
 for j=1:m
if (LM(i)*LM(j))>0 % Si el nodo no está restringido
CONTADOR = CONTADOR+1;
%Este bloque duplica la longitud de los vectores cuando
ya se han llenado los espacios que estaban predefinidos
if CONTADOR>LENGTH_K
% Duplica la longitud del vector K_ROW llenando lo
nuevo con ceros
K_ROW(2*CONTADOR)= 0;
% Duplica la longitud del vector K_COL llenando lo
nuevo con ceros
K_COL(2*CONTADOR)= 0;
% Duplica la longitud del vector K_VAL llenando lo
nuevo con ceros
K_VAL(2*CONTADOR)= 0;
LENGTH_K = 2*CONTADOR;
End
%Fila de la K global donde ira el valor k(i,j)
K_ROW(CONTADOR) = LM(i);
%Columna de la K global donde ira el valor k(i,j)
K_COL(CONTADOR) = LM(j);
% Valor en la K global
K_VAL(CONTADOR) = k(i,j);
end
end
end
end
end

%% ---------------------------------------------------------------------------

%% Creacion de K dispersa y solucion del sistema K*U=F

K = sparse(K_ROW(1:CONTADOR),K_COL(1:CONTADOR),K_VAL(1:CONTADOR

%%Liberar Memoria
clear K_ROW K_COL K_VAL CONTADOR

%%Resolver sistema
U = K\applied_forces; % esto es una function de matlab
R = K*U - applied_forces;
F = R + applied_forces;
disp(F)
%-----------------------------------------------------------------------------

Para encontrar los esfuerzos principales es necesario resolver la siguiente ecuación:

[σ ] = [D ]⋅ [B ]⋅ [u ] (9 .14 )

A continuación se presenta el código en Matlab para el encontrar los esfuerzos

principales
%-----------------------------------------------------------------------------
%% Calculo de esfuerzos principales

%Desplazamientos por nodo

%Se guarda el desplazamiento del nodo respectivo encontrado en U.
displacements_by_node = zeros(n_nodes,n_dims);
for node=1:n_nodes
for dim=1:n_dims
if ID(node,dim)
displacements_by_node(node,dim)=U(ID(node,dim));
end

48
 end
end

%-----------------------------------------------------------------------------

S = zeros(6,n_nodes);
node_count = zeros(n_nodes,1);

%%Tensiones por elemento

s = zeros(6,1);
for element = 1:n_elements
% Se toman los números de los nodos que están presentes en el elemento
nodes = elementsTrian(element,1:Tr);
% Se toman los desplazamientos de los nodos del elemento
a = displacements_by_node(nodes,:)';
% Se hace un vector vertical con los desplazamientos de los nodos
a = a(:);
s(1:3,1) = D{elementsTrian(element,end)} * B{element} * a;
%matriz combined stressed Los eigenvalores se devuelven de menor a mayor
sp = eig([[s(1) s(3)];[s(3) s(2)]]);
vm = sqrt(0.5*((s(1))^2+(s(3))^2+(s(3)-s(1))^2));
s(4:6,1) = [[sp(2) sp(1)]'; vm ];

%suavizado promedio de los esfuerzos por nodo

node_count(nodes) = node_count(nodes) + 1;
for n1 = 1 : Tr
S(:,nodes(n1)) = (S(:,nodes(n1))*(node_count(nodes(n1))-1)
+ s)/node_count(nodes(n1));
end
end
%-----------------------------------------------------------------------------

49
Metodología

En la primera etapa del proyecto se construyo una función en GID que permite por
medio de ficheros almacenar la información de la maya generada por GID, en el primer
fichero se encontrara la cantidad de nodos que posee la maya, cuantos elementos son
triangulares, cuantos cuadriláteros, cuantos tetraedros y por ultimo el tipo de material;
en el segundo fichero se encontrara las restricciones y los momentos para cada nodo
en x, y, z; en los ficheros 3, 4 y 5 se presenta el numero de elemento, números de
nodos de este elemento y los nodos i, j, k, l, m, n, o, p, según sea el caso; en el fichero
sexto se encuentra las coordenadas de los nodos; en el séptimo numero de materiales
y por ultimo en el fichero octavo se encuentra las fuerzas y los momentos en x, y ,z.

Después de esto se construyo una plataforma en Matlab que permite leer y almacenar
la información que se encuentra en cada uno de los ficheros.

Se procedió entonces a realizar la implementación en Matlab de los siguientes

elementos

Elemento triangular de 3 nodos

Elemento triangular de 6 nodos

Elemento cuadrático de 4 nodos

Elemento cuadrático de 9 nodos

Elemento tetraedro de 4 nodos

Elemento hexaédrico de 8 nodos

A partir de esto se construyo una función en Matlab que permite ensamblar la matriz
de rigidez global de una manera mas eficiente (utilizando matrices dispersas), para así
poder encontrar la solución a nuestro problema ( KU = F ).

Trabajo futuro

Después de determinar numéricamente cuales son las tenciones y deformaciones a la

que la pieza es sometida se requiere realizar una interfaz entre Matlab y GID que
permita realizar la apreciación de los resultados de forma visual.

Después que el software este refinado realizarle algunas pruebas y lograr con esto
optimizar el tiempo de computo y o la memoria requerida en su ejecución.

50
Conclusiones

Se exploro los sistemas de almacenamiento de matrices dispersas optimizados

en Matlab.

Se estudiaron varios algoritmos FEA y se programaron los elementos planos y

de volumen.

Se implemento una función que lee la información y determina con esta que
tipo de elementos requiere utilizar.

Se implementó un software inicial que permite calcular problemas hasta de

75000 grados de libertad.

Se probó que los resultados obtenidos coinciden con la solución analítica de un

problema sencillo.

Se debe mejorar el manejo de memoria del programa.

Se logro unificar una estructura en el programa que permite realizar mejoras y

correcciones de una forma fácil, para las personas que intervienen en el
proyecto.

Se desarrolló una aplicación real y se resolvió con el software.

51
Referencias

Pepper, Darrell y Heinrich, Juan. (1992), The Finite Element Method Basic
Concepts and Applications. United States of America: Hemisphere Publishing
Corporation.

Bathe, Klaus Jürgen. (1982), Finite Element Procedures. New Jersey:

Prentice Hall.

Zienkiewicz,O. C. y Leroy , Robert. (2000),The Finite Element Method: Solid

mechanics, Barcelona, España: Printed and bound by MPG Books Ltd.

Fornóns, José-María. (1982), El método de los elementos finitos en la

ingeniería de estructuras. España. Talleres gráficos.

P. G. Ciarlet (1978): The Finite Element Method for Elliptic Problems, North-
Holland, Amsterdam.

P. G. Ciarlet (1991): Basic error estimates for elliptic problems, en Handbook of

Numerical Analysis (Vol II) J.L. Lions y P. G. Ciarlet (ed.), North-Holland,
Amsterdam.

Paginas web

• http://gid.cimne.upc.es/index.html (manual)

• www.mathworks.com

52