1 

Introducción al Método de los 
Elementos Finitos 

1.1. SISTEMAS DISCRETOS Y SISTEMAS CONTINUOS 
Al  efectuar  una  clasificación  de  las  estructuras,  suelen  dividirse  en  discretas  o 
reticulares  y  continuas.  Las  primeras  son  aquéllas  que  están  formadas  por  un 
ensamblaje  de  elementos  claramente  diferenciados  unos  de  otros  y  unidos  en  una 
serie de puntos concretos, de tal manera que el sistema total tiene forma de malla o 
retícula.  La  característica  fundamental  de  las  estructuras  discretas  es  que  su 
deformación  puede  definirse  de  manera  exacta  mediante  un  número  finito  de 
parámetros,  como  por  ejemplo  las  deformaciones  de  los  puntos  de  unión  de  unos 
elementos  y  otros.  De  esta  manera  el  equilibrio  de  toda  la  estructura  puede 
representarse  mediante  las  ecuaciones  de  equilibrio  en  las  direcciones  de  dichas 
deformaciones. 

Figura 1.1  Estructura reticular discreta y estructura continua 

Como  contrapartida,  en  los  sistemas  continuos  no  es  posible  separar,  a  priori,  el 
sistema  en  un  número  finito  de  elementos  estructurales  discretos.  Si  se  toma  una 
parte  cualquiera  del  sistema,  el  número  de  puntos  de  unión  entre  dicha  parte  y  el 
resto de la estructura es infinito, y es por lo tanto imposible utilizar el mismo método 
que en los sistemas discretos, pues los puntos de unión entre los distintos elementos, 
que allí aparecían de manera natural, no existen ahora. 
Las  estructuras  continuas  son  muy  frecuentes  en  ingeniería,  como  por  ejemplo: 
bastidores de máquinas, carrocerías de vehículos, losas de cimentación de edificios, 
vasijas  de  reactores,  elementos  de  máquinas  (bielas,  poleas,  carcasas...),  y  para  su 
análisis  es  necesario  disponer  de  un  método  que  tenga  en  cuenta  su  naturaleza 
continua. 
Hasta la llegada del Método de los Elementos Finitos (MEF), los sistemas continuos se 
abordaban  analíticamente,  pero  por  esa  vía  sólo  es  posible  obtener  solución  para 
sistemas  con  geometría  muy  sencilla,  y/o  con  condiciones  de  contorno  simples. 
También  se  han  utilizado  técnicas  de  diferencias  finitas,  pero  éstas  plantean 
problemas cuando los contornos son complicados. 
Como  precursores  del  MEF  debe  citarse  a  Argyris  y  Kelsey  (Stuttgart,  1955)  y  a 
Turner,  Clough,  Martin  y  Topp  (Boeing, 1956), aunque con posterioridad el número 
de  autores  en  el  campo  del  MEF  ha  sido  enorme,  siendo  uno  de  los  campos  de  la 
ingeniería a los que más esfuerzos de investigación se han dedicado. 

1.2. HIPÓTESIS DE DISCRETIZACIÓN 
En  una  estructura  discreta,  su  deformación  viene  definida  por  un  número  finito  de 
parámetros  (deformaciones  y/o  giros),  que  juntos  conforman  el  vector  de 
deformaciones Δ, y la estructura tiene tantas formas de deformarse como términos 
tenga dicho vector. Un medio continuo tiene infinitas formas posibles de deformarse, 
independientes  unas  de  otras,  ya  que  cada  punto  puede  desplazarse  manteniendo 
fijos cualquier número finito de los puntos restantes, por grande que sea este último. 
Por lo tanto la configuración deformada de la estructura no puede venir dada por un 
vector finito Δ como el anterior, sino que es una función vectorial u, que indica cuáles 
son las deformaciones de cualquier punto, y que tiene tres componentes escalares: 
⎪⎧⎪u(x , y, z )⎪⎫⎪
⎪ ⎪
  u = ⎪⎨ v(x , y, z ) ⎪⎬   (1.1) 
⎪⎪ ⎪
⎪⎪w(x , y, z )⎪⎪⎪
⎩ ⎭
Esta función es la solución de la ecuación diferencial que gobierna el problema, y si 
éste  está  bien  planteado,  cumplirá  las  condiciones  de  contorno  impuestas,  pero  en 
principio  no  puede  asegurarse  que  esta  función  u  tenga  una  expresión  analítica 
manejable,  ni  siquiera  que  pueda  calcularse.  Por  lo  tanto  la  función  u  no  podrá 
conocerse en general. 
Para  resolver  este  problema,  el  Método  de  los  Elementos  Finitos  recurre  a  la 
hipótesis de discretización, que se basa en lo siguiente: 
• El continuo se divide por medio de líneas o superficies imaginarias en una serie de 
regiones  contiguas  y  disjuntas  entre  sí,  de  formas  geométricas  sencillas  y 
normalizadas, llamadas elementos finitos. 
 

• Los  elementos  finitos  se  unen  entre  sí  en  un  número  finito  de  puntos,  llamados 
nudos. 
• Los desplazamientos de los nudos son las incógnitas básicas del problema, y éstos 
determinan unívocamente la configuración deformada de la estructura. Sólo estos 
desplazamientos nodales se consideran independientes. 
• El desplazamiento de un punto cualquiera, viene unívocamente determinado por 
los  desplazamientos  de  los  nudos  del  elemento  al  que  pertenece  el  punto.  Para 
ello se definen para cada elemento, unas funciones de interpolación que permiten 
calcular  el  valor  de  cualquier  desplazamiento  interior  por  interpolación  de  los 
desplazamientos nodales. Estas funciones de interpolación serán de tal naturaleza 
que se garantice la compatibilidad de deformaciones necesaria en los contornos 
de unión entre los elementos. 
• Las  funciones  de  interpolación  y  los  desplazamientos  nodales  definen 
unívocamente  el  estado  de  deformaciones  unitarias  en  el  interior  del  elemento. 
Éstas,  mediante  las  ecuaciones  constitutivas  del  material  definen  el  estado  de 
tensiones en el elemento y por supuesto en sus bordes. 
• Para cada elemento, existe un sistema de fuerzas concentradas en los nudos, que 
equilibran a las tensiones existentes en el contorno del elemento, y a las fuerzas 
exteriores sobre él actuantes. 
Los  dos  aspectos  más  importantes  de  esta  hipótesis,  sobre  los  que  hay  que  hacer 
hincapié son: 
• La  función  solución  del  problema  u  es  aproximada  de  forma  independiente  en 
cada  elemento.  Para  una  estructura  discretizada  en  varios  elementos,  pueden 
utilizarse funciones de interpolación distintas para cada uno de ellos, a juicio del 
analista,  aunque  deben  cumplirse  ciertas  condiciones  de  compatibilidad  en  las 
fronteras entre los elementos. 
• La  función solución es aproximada dentro de cada elemento, apoyándose en un 
número finito (y pequeño) de parámetros, que son los valores de dicha función en 
los nudos que configuran el elemento y a veces sus derivadas. 
Esta hipótesis de discretización es el pilar básico del MEF, por lo que se suele decir de 
éste,  que  es  un  método  discretizante,  de  parámetros  distribuidos.  La  aproximación 
aquí indicada se conoce como la formulación en desplazamiento. 
Claramente  se  han  introducido  algunas  aproximaciones.  En  primer  lugar  no  es 
siempre  fácil  asegurar  que  las  funciones  de  interpolación  elegidas  satisfarán  al 
requerimiento de continuidad de desplazamientos entre elementos adyacentes, por 
lo  que  puede  violarse  la  condición  de  compatibilidad  en  las  fronteras  entre  unos  y 
otros.  En  segundo  lugar  al  concentrar  las  cargas  equivalentes  en  los  nudos,  las 
condiciones  de  equilibrio  se  satisfarán  solamente  en  ellos,  y  no  se  cumplirán 
usualmente en las fronteras entre elementos. 
El proceso de discretización descrito tiene una justificación intuitiva, pero lo que de 
hecho se sugiere es la minimización de la energía potencial total del sistema, para un 
campo  de  deformaciones  definido  por  el  tipo  de  elementos  utilizado  en  la 
discretización. 
Con  independencia  de  que  más  adelante  se  estudien  en  detalle,  se  representan  a 
continuación algunos de los elementos más importantes. 
• Elasticidad unidimensional 

Figura 1.2  Elementos para elasticidad unidimensional 

• Elasticidad bidimensional 

Figura 1.3  Elementos para elasticidad bidimensional 

• Elasticidad tridimensional 

Figura 1.4  Elementos para elasticidad tridimensional 
 

• Elasticidad con simetría de revolución 

Figura 1.5  Elemento axisimétrico 

• Vigas 

Figura 1.6  Elemento viga 

• Flexión de placas planas 

Figura 1.7  Elementos placa plana 

 
• Cáscaras laminares curvas 

Figura 1.8  Elemento cáscara curva 

1.3. FUNCIONES DE INTERPOLACIÓN 
Consideremos  un  elemento  finito  cualquiera,  definido  por  un  número  de  nudos  n. 
Para facilitar la exposición se supondrá un problema de elasticidad plana. Un punto 
cualquiera  del  elemento  tiene  un  desplazamiento  definido  por  un  vector  u,  que  en 
este caso tiene dos componentes: 
⎪⎧u(x , y )⎪⎫⎪
  u = ⎪⎨ ⎬  (1.2) 
⎪⎪v(x , y )⎪⎪
⎩ ⎭
Los nudos del elemento tienen una serie de grados de libertad, que corresponden a 
los valores que adopta en ellos el campo de desplazamientos, y que forman el vector 
denominado  δe . Para el caso plano este vector es: 
T
  δe = ⎡⎢U 1 V1 U 2 V2 ... U n Vn ⎤⎥   (1.3) 
⎣ ⎦
En  este  ejemplo  se  supone  que  como  deformaciones  de los nudos se emplean sólo 
los  desplazamientos,  pero  no  los  giros,  lo  cual  es  suficiente  para  elasticidad  plana, 
como  se  verá  más  adelante.  En  otros  elementos  (p.e.  vigas  o  cáscaras)  se  emplean 
además los giros. 
 

Figura 1.9  Deformaciones en un elemento finito 

El campo de deformaciones en el interior del elemento se aproxima haciendo uso de 
la hipótesis de interpolación de deformaciones: 
  u = ∑ Ni Ui v = ∑ N iVi   (1.4) 

donde  Ni  son  las  funciones  de  interpolación  del  elemento,  que  son  en  general 
funciones de las coordenadas x,y. Nótese que se emplean las mismas funciones para 
interpolar los desplazamientos u y v, y que ambos desplazamientos se interpolan por 
separado, el campo u mediante las Ui y el campo v mediante las Vi. Es decir que la 
misma  Ni define  la  influencia  del  desplazamiento  del  nudo  i  en  el  desplazamiento 
total del punto P, para las dos direcciones x e y. 
La interpolación de deformaciones (1.4) puede ponerse en la forma matricial general: 
  u = N δe   (1.5) 

La  matriz  de  funciones  de  interpolación  N  tiene  tantas  filas  como  desplazamientos 
tenga  el  punto  P  y  tantas  columnas  como  grados  de  libertad  haya  entre  todos  los 
nudos del elemento.  
Las  funciones  de  interpolación  son  habitualmente  polinomios,  que  deben  poderse 
definir  empleando  las  deformaciones  nodales  del  elemento.  Por  lo  tanto  se  podrán 
usar polinomios con tantos términos como grados de libertad tenga el elemento. 
Para problemas de elasticidad la estructura de esta matriz es normalmente del tipo: 
⎡N 0 N2 0 ... 0 N n 0 ⎤⎥
  N = ⎢⎢ 1   (1.6) 
⎢⎣ 0 N 1 0 N2 0 ... 0 N n ⎥⎥
⎦
Sin  embargo,  el  aspecto  de  esta  matriz  puede  ser  distinto  para  otros  elementos, 
como las vigas o las placas a flexión. 
Las  funciones  de  interpolación  están  definidas  únicamente  para  el  elemento,  y  son 
nulas  en  el  exterior  de  dicho  elemento.  Estas  funciones  tienen  que  cumplir 
determinadas  condiciones  y  aunque  éstas  se  verán  en  detalle  más  adelante,  con  la 
expresión anterior se puede deducir que la función de interpolación Ni debe valer 1 
en el nudo i y 0 en los restantes nudos. Esta condición resulta evidente si se tiene en 
cuenta  que  los  términos  del  vector  δe   son  grados  de  libertad  y  por  lo  tanto  son 
independientes, y deben poder adoptar cualquier valor. 

Figura 1.10  Función de interpolación 

1.4. CRITERIOS DE CONVERGENCIA 
Antes de estudiar los criterios para garantizar la convergencia en el MEF es necesario 
definir dicho concepto, en el ámbito del MEF. Se dice que un análisis por el MEF es 
convergente si al disminuir el tamaño de los elementos, y por lo tanto aumentar el 
número  de  nudos  y  de  elementos,  la  solución  obtenida  tiende  hacia  la  solución 
exacta. 
Hay que indicar que en el análisis por el MEF, se introducen, además de la hipótesis 
de  discretización,  otras  aproximaciones,  que  son  fuentes  de  error  en  la  solución: 
integración  numérica,  errores  de  redondeo  por  aritmética  finita...  El  concepto  de 
convergencia  aquí  analizado  se  refiere  solamente  a  la  hipótesis  de  discretización, 
prescindiendo  de  los  otros  errores,  que  deben  ser  estudiados  aparte,  y  cuyo  valor 
debe en todo caso acotarse. 
Las funciones de interpolación elegidas para representar el estado de deformación de 
un medio continuo deben satisfacer una serie de condiciones, a fin de que la solución 
obtenida por el MEF, converja hacia la solución real. 
Criterio 1 
Las funciones de interpolación deben ser tales que cuando los desplazamientos de los 
nudos  del  elemento  correspondan  a  un  movimiento  de  sólido  rígido,  no  aparezcan 
tensiones en el elemento. 
Este  criterio  se  puede  enunciar  también  de  forma  más  sencilla:  las  funciones  de 
interpolación  deben  ser  capaces  de  representar  los  desplazamientos  como  sólido 
rígido, sin producir tensiones en el elemento. 
Esta condición es evidente, pues todo sólido que se desplaza como un sólido rígido, 
no sufre ninguna deformación ni por lo tanto tensión. Sin embargo adoptando unas 
 

funciones de interpolación incorrectas, pueden originarse tensiones al moverse como 
sólido rígido. 
Por ejemplo en la figura 1.11, los elementos del extremo se desplazan como un sólido 
rígido, al no existir tensiones más allá de la fuerza aplicada. 

Figura 1.11  Deformación de sólido rígido 

Empleando  la  formulación  desarrollada  más  adelante,  si  se  aplican  unas 
deformaciones  en  los  nudos  de  valor  δR  que  representan  un  movimiento  de  sólido 
rígido, las deformaciones unitarias en el interior del elemento son: 
  εR = B δ R   (1.7) 

Según  este  criterio,  las  tensiones  correspondientes  deben  ser  nulas  en  todo  punto 
del elemento: 
  σ = D εR = D B δ R = 0   (1.8) 

Criterio 2 
Las funciones de interpolación deben ser tales que cuando los desplazamientos de los 
nudos  correspondan  a  un  estado  de  tensión  constante,  este  estado  tensional  se 
alcance en realidad en el elemento. 
Claramente,  a  medida  que  los  elementos  se  hacen  más  pequeños,  el  estado  de 
tensiones que hay en ellos se acerca al estado uniforme de tensiones. Este criterio lo 
que exige es que los elementos sean capaces de representar dicho estado de tensión 
constante. 
Se observa que este criterio de hecho es un caso particular del criterio 1, ya que un 
movimiento como sólido rígido (con tensión nula) es un caso particular de un estado 
de  tensión  constante.  En  muchas  ocasiones  ambos  criterios  se  formulan  como  un 
sólo criterio.  
A los elementos que satisfacen los criterios 1 y 2 se les llama elementos completos. 
Criterio 3 
Las funciones de interpolación deben ser tales que las deformaciones unitarias que se 
produzcan  en  las  uniones  entre  elementos  deben  ser  finitas.  Esto  es  lo  mismo  que 
decir que debe existir continuidad de desplazamientos en la unión entre elementos 
aunque puede haber discontinuidad en las deformaciones unitarias (y por lo tanto en 
las tensiones, que son proporcionales a ellas).  
La figura 1.12 ilustra las posibles situaciones, para un caso unidimensional donde la 
única  incógnita  es  el  desplazamiento  u  en  la  dirección  x.  En  la  situación  de  la 
izquierda  existe  una  discontinuidad  en  el  desplazamiento  u,  que  da  lugar  a  una 
deformación unitaria infinita: esta situación no está permitida por el criterio 3. En la 
situación mostrada a la derecha la deformación es continua, aunque la deformación 
unitaria no lo es: esta situación está permitida por el criterio 3. 

Figura 1.12  Criterio de convergencia 3 en una dimensión 

Este criterio debe cumplirse para poder calcular la energía elástica U almacenada en 
toda la estructura, como suma de la energía de todos los elementos.  
  U = 1
2 ∫σ
T
ε dv = 1
2 ∑∫ σ T
ε dv + U cont   (1.9) 
V e ve

donde  el  sumando  Ucont  representa  la  energía  elástica  acumulada  en  el  contorno 
entre los elementos, que debe ser nula. 
 

Si  este  requerimiento  no  se  cumple,  las  deformaciones  no  son  continuas  y  existen 
deformaciones  unitarias  ε  infinitas  en  el  borde  entre  elementos.  Si  la  deformación 
unitaria en el contorno es infinita, la energía que se acumula en él es: 
  U cont = 1
2 ∫ σT (∞) dv = indeterminado   (1.10) 
v =0

ya que, aunque el volumen de integración (volumen del contorno) es nulo, su integral 
puede  ser  distinta  de  cero,  con  lo  que  se  almacena  energía  en  los  bordes  entre 
elementos, lo cual no es correcto. 
Sin  embargo,  si  la  deformación  unitaria  en  el  contorno  es  finita  (aunque  no  sea 
continua entre los elementos unidos), la energía que se acumula es: 
  U cont = 1
2 ∫ σT εcont dv = 0   (1.11) 
v =0

En  el  caso  plano  o  espacial  este  requerimiento  aplica  a  la  componente  del 
desplazamiento perpendicular al borde entre elementos, ya que ésta es la única que 
acumula energía.  

Figura 1.13  Compatibilidad de desplazamientos en el contorno 

Este criterio puede expresarse de manera más general diciendo que en los contornos 
de los elementos deben ser continuas las funciones de interpolación y sus derivadas 
hasta un orden n-1, siendo n el orden de las derivadas existentes en la expresión de 
la energía potencial Π del sistema. Es decir que las funciones de interpolación deben 
ser continuas de tipo Cn-1. 
El orden n de las derivadas existentes en la energía potencial Π del sistema, siempre 
es  la  mitad  del  orden  de  la  ecuación  diferencial  que  gobierna  el  problema  m 
(n=m/2). 
Para elementos de tipo celosía o viga, este requerimiento es fácil de cumplir pues la 
unión  entre  elementos  se  hace  en  puntos  discretos,  y  se  usan  los  mismos 
desplazamientos y giros para todos los elementos que se unen en un nudo. 
Para  elasticidad  plana  la  ecuación  diferencial  es  de  orden  m=2,  con  lo  que  energía 
potencial es de orden n=1. En efecto esta última se expresa en términos de ε, que 
son  las  derivadas  primeras  de  las  deformaciones.  Luego  las  funciones  de 
interpolación deben ser continuas en los contornos de tipo C0, es decir no se exige 
continuidad a la derivada de la función de interpolación en el contorno del elemento, 
sino sólo a la propia función. 
Para problemas de flexión de vigas y de placas delgadas, la ecuación diferencial es de 
orden m=4, luego la energía potencial es de orden n=2. Por lo tanto las funciones de 
interpolación elegidas deben ser continuas C1 en el contorno del elemento: tanto la 
función como su derivada primera deben ser continuas. 
A los elementos que cumplen este tercer criterio se les llama compatibles. 

Resumen de los tres criterios

Los  criterios  1  y  2  pueden  agruparse  de  manera  más  matemática  diciendo  que  las 
funciones de interpolación deben permitir representar cualquier valor constante de 
su  derivada  n‐sima  en  el  interior  del  elemento  (siendo  n  el  orden  de  la  derivada 
necesaria  para  pasar  de  las  deformaciones  a  las  deformaciones  unitarias,  que  es  el 
orden de derivación de las deformaciones en el potencial). 
Esto puede comprobarse mediante el siguiente sencillo razonamiento: los criterios 1 
y  2  obligan  a  representar  cualquier  valor  constante  (incluido  el  valor  nulo)  de  la 
tensión σ, lo cual equivale a representar cualquier valor constante (incluido el valor 
nulo) de la deformación unitaria ε. Pero la deformación unitaria es la derivada n‐sima 
de la deformación, luego en consecuencia es necesario poder representar cualquier 
valor  constante  de  dicha  derivada.  Esto  se  satisface  siempre,  si  se  adoptan  como 
funciones de interpolación, polinomios completos de orden n como mínimo. 
Por ejemplo, si se adopta una función lineal N=Ax+B, sólo es válida para problemas 
de elasticidad (n=1), ya que se representa cualquier valor constante de la derivada n‐
sima dN/dx=A. Sin embargo, no vale para problemas de flexión de vigas ni de placas 
delgadas (n=2), pues siempre es d2N/dx2=0, es decir que no se puede representar 
cualquier valor constante de la derivada segunda. 
El criterio 3 exige que las deformaciones unitarias sean finitas en el contorno entre 
los elementos. Como estas deformaciones unitarias son las derivadas n‐simas de las 
deformaciones, lo que se exige es que haya continuidad de las deformaciones y sus 
derivadas  hasta  orden  n-1  en  el  contorno  del  elemento.  Esto  es  equivalente  a 
imponer la compatibilidad de desplazamientos en el contorno. 
Como resumen de los tres criterios, para problemas de elasticidad (n=1) es necesario 
emplear  polinomios  completos  de  orden  1,  con  continuidad  C0  entre  ellos  para 
garantizar  la  convergencia.  Es  suficiente  con  usar  funciones  del  tipo  lineal,  que 
aproximan  la  solución  mediante  una  línea  quebrada,  aunque  se  produzcan 
discontinuidades  en  las  tensiones  entre  los  elementos,  como  se  indica  en  la  figura 
1.14. 
 

Figura 1.14  Compatibilidad en elasticidad unidimensional 

Para  problemas  de  flexión  de  vigas  y  placas,  (n=2)  es  necesario  emplear  como 
mínimo polinomios de grado 2, con continuidad C1 entre ellos, es decir que hay que 
garantizar  la  continuidad  de  la  flecha  y  el  giro  entre  los  elementos.  En  la  práctica, 
para la flexión de vigas planas se usan 4 parámetros para ajustar la solución (flecha y 
giro en cada extremo) por lo que el tipo de funciones empleadas son polinomios de 
grado 3. 

P
1 v 2 3

v
2

V2 V3
V1

Figura 1.15  Compatibilidad en flexión de vigas 

Con  respecto  a  la  velocidad  de  convergencia  se  pueden  resumir  las  siguientes 
conclusiones de los análisis efectuados. Si se utiliza una discretización uniforme con 
elementos  de  tamaño  nominal  h,  y  se  usa  para  interpolar  los  desplazamientos  un 
polinomio  completo  de  grado  c  (que  representa  exactamente  variaciones  del 
desplazamiento de dicho grado), el error local en los desplazamientos se estima que 
es del orden 0(hc+1). 
Respecto a las tensiones, son las derivadas n‐simas de los desplazamientos, luego el 
error en ellas es de orden 0(hc+1-m). 
 
 2
Ecuaciones generales

2.1. CAMPO DE DEFORMACIONES

El campo de deformaciones en un punto cualquiera del dominio está definido por un
vector u que tiene tantas componentes como deformaciones existen en el dominio.
Para el caso de un problema espacial es:

u(x , y, z )
u v(x , y, z ) (2.1)
w(x , y, z )

Si se considera un elemento finito cualquiera, el campo de deformaciones en su

interior se aproxima, haciendo uso de la hipótesis de interpolación, como un
promedio ponderado de las deformaciones en cada uno de los n nudos del elemento,
siendo los factores de ponderación las funciones de interpolación:
u Ni Ui v NV
i i
w NW
i i
(2.2)

Esta interpolación puede ponerse en forma matricial:

e
u N (2.3)
e
donde es el vector de todas las deformaciones nodales del elemento (figura 2.1):
T
e
U 1 V1 W1 U 2 V2 W2 ... U n Vn Wn (2.4)
 W1
w
v
V1
U1
u

Figura 2.1 Deformaciones en un elemento finito

La matriz de funciones de interpolación N tiene tres filas y tantas columnas como

grados de libertad haya entre todos los nudos del elemento. La estructura de esta
matriz siempre es del tipo:

N1 0 0 ... ... N n 0 0
N 0 N1 0 ... ... 0 Nn 0 (2.5)
0 0 N 1 ... ... 0 0 Nn

2.2. DEFORMACIONES UNITARIAS

Las deformaciones unitarias en un punto cualquiera del elemento, con la suposición
de pequeñas deformaciones, son:
u
x
v
x
y
y w
z z
(2.6)
xy u v
yz y x
zx
v w
z y
w u
x z
Se pueden poner en la forma matricial siguiente:
 0 0
x
0 0
x
y
y
0 0 u
z z
v u (2.7)
xy
0 w
yz y x
zx
0
z y

0
z x
En esta expresión se identifica el operador matricial  que permite pasar de las
deformaciones de un punto u a las deformaciones unitarias . Este operador tiene
tantas filas como deformaciones unitarias haya en el problema y tantas columnas
como componentes tenga el campo de desplazamientos u.
Sustituyendo las deformaciones u en función de las deformaciones nodales, mediante
las funciones de interpolación, se obtiene:
e
u N (2.8)

En esta relación se identifica la matriz B:

B N (2.9)

tal que se cumple que:

e
B (2.10)

Esta matriz B relaciona las deformaciones de los nudos del elemento e con las
deformaciones unitarias en un punto interior cualquiera del elemento. Por lo tanto B
representa el campo de deformaciones unitarias que se supone existe en el interior
del elemento finito, como consecuencia de la hipótesis de interpolación de
deformaciones efectuada, y juega un papel fundamental en el método de los
elementos finitos.
Dada la estructura de la matriz N, la matriz B se puede poner siempre en la forma:

N1 0 0 ... ... N n 0 0
B N 0 N1 0 ... ... 0 Nn 0 (2.11)
0 0 N 1 ... ... 0 0 Nn

B B1 B2 ... Bn (2.12)
Cada una de las matrices Bi tiene la forma siguiente:
Ni
0 0
x
Ni
0 0
y
Ni 0 0 Ni
0 0
z
Bi 0 Ni 0 (2.13)
Ni Ni
0 0 Ni 0
y x
Ni Ni
0
z y
Ni Ni
0
z x
Aunque el valor de B se ha obtenido para el caso de elasticidad tridimensional, su
valor en función de  y N es totalmente general para otros tipos de problemas de
elasticidad, como flexión de placas, problemas de revolución, etc.

2.3. ESTADO DE TENSIONES. ECUACIÓN CONSTITUTIVA

Las tensiones en un punto cualquiera del dominio están definidas por el tensor de
tensiones en dicho punto, cuya expresión general es:
x

z
(2.14)
xy

yz

zx

Asimismo se conoce la ecuación constitutiva del material que forma el dominio, y que
relaciona las tensiones con las deformaciones unitarias. Para un material elástico
lineal esta ecuación constitutiva se puede poner en la forma:
D( 0
) 0
(2.15)

Siendo:
 D la matriz elástica, que para un material elástico lineal es constante y depende
de sólo dos parámetros: el módulo de elasticidad E y el módulo de Poisson .
 0 el vector de las deformaciones unitarias iniciales existentes en el material en el
punto considerado, que deben ser conocidas. Las más habituales son las debidas
 a las temperaturas, aunque pueden incluirse en ellas las debidas a los errores de
forma, etc.
 0 las tensiones iniciales presentes en el material, que normalmente son tensiones
residuales debidas a procesos anteriores sobre el material (p.e. tratamiento
térmico) y que por lo tanto son conocidas.
Las expresiones particulares de la matriz elástica D y de los vectores 0 y 0 dependen
del tipo de problema considerado y serán estudiadas en cada caso particular.

2.4. ECUACIÓN DE EQUILIBRIO DE UN ELEMENTO

Una vez que han quedado establecidas las expresiones que relacionan los
desplazamientos, las deformaciones unitarias y las tensiones, en función de los
desplazamientos de los nudos, se está ya en condiciones de calcular las ecuaciones de
equilibrio de un elemento finito. Si se considera un elemento finito cualquiera, las
fuerzas que actúan sobre él, en el caso más general, son las siguientes (figura 2.2):
 Fuerzas exteriores de volumen aplicadas en el interior del elemento qv, que son en
general variables dentro del elemento, y tienen tantas componentes como
desplazamientos haya en cada punto.
 Fuerzas exteriores de superficie aplicadas en el contorno libre del elemento qs,
que son en general variables a lo largo del contorno, y tienen tantas componentes
como desplazamientos tenga cada punto del contorno. Al contorno sobre el que
actúan las fuerzas de superficie se le denomina s.
 Fuerzas interiores qc, aplicadas en la superficie del contorno de unión del
elemento con los elementos vecinos, que son desconocidas. A dicho contorno de
unión se le denomina c.
 Fuerzas exteriores puntuales aplicadas sobre los nudos del elemento PNe .

PN

qs
qc qv

Figura 2.2 Fuerzas actuantes sobre un elemento finito

El trabajo virtual que producen estas fuerzas es:

We uT qv dv uT qs ds uT qc ds eT
PNe (2.16)
v s c

Donde u es una variación virtual del campo de deformaciones u y e es la variación

correspondiente a los grados de libertad de los nudos. Durante estas variaciones, las
fuerzas exteriores se mantienen constantes.
Aplicando el principio de los trabajos virtuales se obtiene que para que haya
equilibrio, el trabajo virtual de las fuerzas debe ser igual a la variación de la energía
elástica U acumulada en el elemento, para cualquier u:

We T
dv Ue (2.17)
v

Donde  es la variación en las deformaciones unitarias producida por la variación en

las deformaciones u. Por lo tanto la ecuación de equilibrio del elemento es:

uT qv dv uT qs ds uT qc ds eT
PNe T
dv (2.18)
v s c v

Aplicando la hipótesis de interpolación de deformaciones, la variación virtual del

campo de deformaciones es:
e
u N (2.19)

La variación de las deformaciones unitarias se relaciona con la variación de las

deformaciones nodales a través de la matriz B:
e
B (2.20)

Sustituyendo las variaciones u y  en la expresión (2.18) se obtiene la ecuación de

equilibrio aproximada mediante la hipótesis de interpolación de deformaciones:

eT
NT q v dv NT qs ds NT qc ds PNe eT
BT dv (2.21)
v s c v

Considerando que esta ecuación se debe cumplir para cualquier variación arbitraria
de las deformaciones, se obtiene:

NT qv dv NT qs ds NT qc ds PNe BT dv (2.22)
v s c v

Esta ecuación representa el equilibrio del elemento finito. Antes de seguir

desarrollándola, la integral debida a las fuerzas distribuidas qc sobre el contorno de
unión (desconocidas) se sustituye por:

NT qc ds Pce (2.23)
c
 PC
qc

Figura 2.3 Fuerzas de conexión entre elementos

Donde Pce son unas fuerzas que están aplicadas sobre los nudos del elemento, y que
son equivalentes a las fuerzas distribuidas aplicadas sobre los contornos de unión con
los elementos vecinos. Ambas fuerzas producen el mismo trabajo virtual. La ecuación
de equilibrio del elemento queda finalmente:

NT qv dv NT qs ds Pce PNe BT dv (2.24)

v s v

Sustituyendo en ella el valor de la tensión mediante la ecuación constitutiva (2.15) se

obtiene:

NT qv dv NT qs ds Pce PNe BT (D D 0 0 )dv (2.25)

v s v

Sustituyendo a continuación el valor de la deformación unitaria en función de la

matriz B se obtiene:

NT qv dv NT qs ds Pce PNe BT (D B e
D 0 0 )dv (2.26)
v s v

Reordenando los distintos términos se llega a:

BT D B dv e

v
(2.27)
T T T T e e
N q v dv N qs ds B D 0 dv B 0 dv Pc P N
v s v v

Esta es la ecuación final de equilibrio del elemento finito considerado. En ella se

identifican los siguientes términos:
 Matriz de rigidez del elemento finito. Se trata de una matriz cuadrada simétrica de
tamaño igual al número de grados de libertad del elemento.

Ke BT D B dv (2.28)
v
 Vector de fuerzas nodales equivalentes debido a las fuerzas actuantes por unidad
de volumen (figura 2.4).

Pve NT q v dv (2.29)
v

Pv

qv

Figura 2.4 Fuerzas nodales equivalentes a las fuerzas de volumen

 Vector de fuerzas nodales equivalentes a las fuerzas exteriores de superficie.

Pse NT qs ds (2.30)
s

 Vector de fuerzas nodales equivalentes producidas por las deformaciones iniciales

existentes en el material:

PTe BT D 0
dv (2.31)
v

 Vector de fuerzas nodales equivalentes debidas a las tensiones iniciales existentes

en el material:

Pbe BT 0
dv (2.32)
v

La ecuación de equilibrio del elemento puede ponerse en forma compacta como:

Ke e
Pve Pse PTe Pbe Pce PNe (2.33)

Esta ecuación de equilibrio está referida al sistema de ejes en el que se hayan

definido las coordenadas y las deformaciones de los nudos, y al que lógicamente
también se habrán referido las distintas fuerzas actuantes. En ella son conocidos
todos los términos de carga salvo el debido a las fuerzas distribuidas interiores Pce
que se producen en el contorno de unión con los elementos vecinos.
2.5. ECUACIÓN DE EQUILIBRIO DEL CONJUNTO
La ecuación de equilibrio obtenida para un elemento puede aplicarse a todos y cada
uno de los elementos en que se ha dividido el sistema continuo. De esta manera se
garantiza el equilibrio de todos y cada uno de ellos individualmente, apareciendo en
dichas ecuaciones las fuerzas de conexión entre unos y otros elementos.
Para obtener la ecuación de equilibrio de toda la estructura es necesario además
imponer el equilibrio de las fronteras de unión entre los elementos. En estas
fronteras se han introducido las fuerzas de conexión entre los elementos qc, que a su
vez han dado lugar a las fuerzas nodales equivalentes correspondientes Pc, y que
como se ha visto son energéticamente equivalentes a ellas. Por lo tanto el considerar
el equilibrio de las fronteras es equivalente a considerar el equilibrio de los nudos de
unión entre los elementos.
Si se plantean conjuntamente las ecuaciones de equilibrio de todos los nudos de
unión entre todos los elementos, se obtiene un conjunto de ecuaciones que
representa el equilibrio de toda la estructura. Estas ecuaciones se obtienen por
ensamblado de las ecuaciones de equilibrio de los distintos elementos finitos que la
forman, en la forma:

e
Ke e
e
Pve Pse PTe Pbe PNe e
Pce (2.34)

Donde se ha empleado el símbolo e para indicar el ensamblado de las distintas

magnitudes según los grados de libertad de la estructura.
En este proceso de ensamblado se cancelan todas las fuerzas de conexión entre unos
elementos y los vecinos, pues se trata de fuerzas iguales y de signo contrario:

e
Pce 0 (2.35)

Al ser el sistema lineal, el término de la izquierda puede ponerse siempre como:

e
Ke e
K (2.36)

Siendo:
 K es la matriz de rigidez de la estructura completa, obtenida por ensamblaje de
las matrices de los distintos elementos según los grados de libertad
correspondientes a cada uno.
  es el vector de grados de libertad de toda la estructura, que agrupa a los grados
de libertad de todos los nudos.
De esta manera, la ecuación de equilibrio del conjunto de la estructura queda:
K Pv Ps PT Pb PN (2.37)
En esta ecuación:
 PN es el vector de fuerzas exteriores actuantes sobre los nudos de unión de los
elementos.
 Pv, Ps, PT, Pb, son los vectores de fuerzas nodales equivalentes producidos por
las fuerzas de volumen, de superficie, las deformaciones iniciales y las tensiones
iniciales. Son todos conocidos y se obtienen por ensamblado de los
correspondientes a los distintos elementos, según los nudos y direcciones
adecuados. Respecto al vector Pb hay que decir que si el estado de tensiones
iniciales actuantes sobre la estructura está auto-equilibrado, este vector es nulo.
Esto ocurre normalmente con las tensiones residuales. Sin embargo estas
tensiones no están equilibradas si la estructura se obtiene partir de un material
con unas tensiones auto-equilibradas y se elimina parte de ese material.

2.6. CONDICIONES DE LIGADURA

La ecuación de equilibrio deducida antes representa el equilibrio del conjunto de la
estructura, considerando todos los elementos que la forman, así como todas las
fuerzas exteriores. Para poderla resolver es necesario imponer las condiciones de
ligadura, que indican cómo está sustentada la estructura.
La introducción de estas condiciones de ligadura se efectúa exactamente igual que en
el método de rigidez para estructuras reticulares.

2.7. ENERGÍA POTENCIAL TOTAL

La densidad de energía elástica acumulada en un punto del elemento es:

U 0e T
d (D D 0 0
)T d 1
2
T
D T
0
D T
0
(2.38)
0 0

El potencial total acumulado en un elemento finito cualquiera es igual a la suma de la

energía elástica acumulada en él más el potencial de las fuerzas exteriores V:
e
Ue Ve U 0e dv Ve (2.39)
v

e T T T
1
2 D dv D 0 dv 0 dv
v v v
(2.40)
uT q v dv uT qs ds uT qc ds eT
PNe
v s c

Sustituyendo las deformaciones unitarias  y los desplazamientos u en función de las

deformaciones nodales mediante las matrices B y N se obtiene:
 e eT
1
2 BT D B dv e eT
BT D 0
dv eT
BT 0
dv
v v v

eT
NT qv dv eT
NT qs ds eT
NT qc ds eT
PNe (2.41)
v s c

En esta expresión se identifican la matriz de rigidez del elemento, así como los
distintos vectores de fuerzas nodales equivalentes, con lo que se puede poner en
forma más compacta:
e eT
1
2 Ke e eT
PTe eT
Pbe eT
Pve eT
Pse eT
Pce eT
PNe (2.42)

El potencial total para el medio continuo se obtiene sumando el potencial de cada

uno de los elementos que lo forman:
e
(2.43)
e

Al sumar el potencial de los distintos elementos, se pueden ir ensamblando la matriz

de rigidez y los vectores de fuerzas de los elementos según los distintos grados de
libertad de la estructura, para obtener la siguiente expresión:
T T T T T T
1
2 K PT Pb Pv Ps c PN (2.44)

En ella aparecen la matriz de rigidez de toda la estructura, y los correspondientes

vectores de grados de libertad y fuerzas nodales equivalentes.
El término c corresponde al potencial acumulado en las fronteras de conexión entre
los elementos por las fuerzas de conexión qc. Este término es nulo si las funciones de
interpolación se eligen con la condición de que no se acumule energía en las
fronteras, como así se hace (véase el tercer criterio de convergencia).
El equilibrio de la estructura implica que el potencial total sea estacionario, para
cualquier variación de las deformaciones nodales:

0 0 (2.45)

Se obtiene así la ecuación de equilibrio de toda la estructura:

K PT Pb Pv Ps PN (2.46)
 3
Elasticidad unidimensional

3.1. INTRODUCCIÓN
En el problema unidimensional el dominio continuo que se analiza se extiende según
una única dimensión x, teniendo el material un área variable con dicha coordenada
A(x) (figura 3.1). Como fuerzas exteriores pueden actuar:
 Fuerzas por unidad de volumen q, en sentido axial.
 Fuerzas puntuales aplicadas Fci.
El campo de deformaciones es una función u(x) de una sola variable, que define la
deformación axial de un punto cualquiera del dominio. La deformación unitaria tiene
sólo la componente axial du / dx .

u(x) q x

Figura 3.1 Problema de elasticidad unidimensional

El problema de elasticidad unidimensional no es de gran aplicación práctica, pero su

estudio tiene interés pues permite conocer las peculiaridades del MEF con ejemplos
muy sencillos.
La ecuación diferencial que gobierna el problema se obtiene aplicando el equilibrio
de fuerzas a un elemento diferencial (figura 3.2):
F (F dF ) q Adx (3.1)

dF
qA 0 (3.2)
dx
Sustituyendo el valor de la fuerza en función de la tensión F A , y ésta en función
de la deformación unitaria E se obtiene:
d du
EA qA 0 (3.3)
dx dx

q A dx
F F+dF

dx

Figura 3.2 Equilibrio de fuerzas en elasticidad unidimensional

Las condiciones de contorno pueden ser de dos tipos:

 Desplazamiento conocido en algunos puntos u=uc
 Fuerza aplicada conocida en algunos puntos F=Fci  EA du/dx=Fci
La energía elástica acumulada en el material es:
2
2 du
U 1
2 dv 1
2 E dv 1
2 EA dx (3.4)
dx
El potencial total es:
2
du
U V 1
2 EA dx q u Adx Fci uci (3.5)
dx
El orden de derivación del desplazamiento en el potencial es 1, por lo que se
requieren funciones de interpolación de tipo C0: sólo se exige continuidad a la función
para garantizar la convergencia del método en este problema.

3.2. ELEMENTO DE DOS NUDOS

El elemento más simple para este problema tiene dos nudos, con desplazamientos U1
y U2 en ellos. La interpolación del desplazamiento dentro del elemento es:
u N1 U1 N2 U2 (3.6)

U1
e
u N1 N 2 N (3.7)
U2
 U1 U2

Figura 3.3 Elemento de dos nudos

Con dos parámetros nodales, se puede interpolar una línea recta: u ax b

Particularizando en los nudos 1 y 2 se obtiene:
U1 a x1 b
(3.8)
U2 a x2 b

De estas expresiones se obtienen a y b,

U2 U1 x2 U1 x1U 2
a b (3.9)
L L
Sustituyendo en la expresión inicial y reordenando, se obtiene la ley de interpolación:
x2 x x x1
u U1 U (3.10)
x 2 x1 x 2 x1 2
Las funciones de interpolación son líneas rectas que valen 1 en el nudo i y 0 en el otro
nudo:
x2 x x x1
N1 N2 (3.11)
x 2 x1 x 2 x1

U2
U1
N1 N2
+1 +1

Figura 3.4 Funciones de interpolación. Elemento de dos nudos

La deformación unitaria es:

du e dN 1 dN 2 U 1
N (3.12)
dx dx dx U 2

Se define la matriz B como:

dN 1 dN 2 1 1
B (3.13)
dx dx L L

Siendo L=x2-x1 la longitud del elemento.

La matriz elástica es sencillamente D = [ E ]
La matriz de rigidez resulta ser:
1
L 1 1
K BT D B dv E Adx (3.14)
1 L L
L
Suponiendo E y A constantes se obtiene:
EA 1 1 EA 1 1
K dx (3.15)
L2 1 1 L 1 1

Esta es la matriz de rigidez del elemento de celosía, ya que este elemento finito de
dos nudos coincide con dicho elemento estructural.

3.2.1 Elemento con área variable

Los grados de libertad, las funciones de interpolación y la matriz B son las mismas que
en el elemento de área constante.

U1 A(x) U2

Figura 3.5 Elemento de área variable

La matriz de rigidez es:

1
L 1 1 E 1 1
K E Adx Adx (3.16)
1 L L L2 1 1
L
En la integral se identifica el área media del elemento Am, con lo que se obtiene:
E Am 1 1
K (3.17)
L 1 1

Se obtiene la misma expresión que para el elemento de área constante, pero usando
el área media del elemento.

3.2.2 Tensiones
La tensión en el elemento (no incluyendo el efecto de las temperaturas) es:
 du e 1 1 U1
E E EB E (3.18)
dx L L U2

E
U U1 (3.19)
L 2
Se observa que el elemento produce un campo de tensión constante en todo su
interior.
Además la tensión también es constante si el elemento es de área variable, en
contradicción con la estática, pues lo que debe ser constante en este caso es la fuerza
axial N, y no la tensión. Esto es debido a la hipótesis efectuada de variación lineal de
la deformación u, que no es correcta para un elemento de área variable.

3.3. FUNCIONES DE INTERPOLACIÓN

Se trata de definir las funciones de interpolación para elementos de tipo general, con
dos o más nudos. Para ello resulta muy práctico el utilizar una coordenada  local al
elemento, de tal forma que ésta siempre varíe de -1 en el nudo inicial a +1 en el nudo
final. En este sistema de coordenadas local, el elemento siempre es de longitud 2.

=-1 =0 =+1

Figura 3.6 Coordenadas locales

 Elemento de dos nudos.

Las funciones de interpolación son:
1 1
N1 N1 (3.20)
2 2

1 N1 N2 1

Figura 3.7 Elemento de dos nudos. Funciones de interpolación

 Elemento de tres nudos.

La función N2 es una parábola que debe valer cero en las dos esquinas y 1 en el nudo
central 2, luego su ecuación debe ser del tipo:
 N2 C (1 )(1 ) (3.21)

Pero en =0 debe ser N2(0)=1, de donde se deduce que el coeficiente C debe valer 1.
Por consiguiente la función es:
2
N2 (1 )(1 ) 1 (3.22)

N2

1

1 2 3

Figura 3.8 Elemento de 3 nudos. Función de interpolación del nudo central

La función N1 es una parábola que debe valer cero en los nudos 2 y 3, y debe valer la
unidad en el nudo 1 (figura 3.9), luego debe ser del tipo:
N1 C (1 ) (3.23)

Pero en =-1 debe ser N1(-1)=1, luego el coeficiente C debe ser –1/2. Por
consiguiente la función es:
( 1)
N1 (3.24)
2

N1

1

1 2 3

Figura 3.9 Elemento de 3 nudos. Función de interpolación de nudo esquina

Análogamente, para el nudo 3 se obtiene (figura 3.10):

( 1)
N3 (3.25)
2
Por lo tanto este elemento tiene una ley parabólica para al interpolación del
desplazamiento (figura 3.11).
 N3 U3
+1
U1 U2

1 2 3 1 2 3

Figura 3.10 Elemento de 3 nudos. Figura 3.11 Elemento de 3 nudos.

Interpolación de nudo esquina Campo de deformaciones

 Elemento de cuatro nudos.

Sus funciones de interpolación se obtienen de forma análoga al elemento de tres
nudos.
1
N1 (1 )(1 3 )(1 3 ) (3.26)
16
9
N2 (1 )(1 )(1 3 ) (3.27)
16
La obtención de las funciones de interpolación siguiendo este método requiere en
algunas ocasiones cierta dosis de buena suerte, y fue bautizada por Zienkiewicz como
formulación “serendipity”, inspirándose en la obra de Horace Walpole “Los Príncipes
de Serendip” (1754).

1 2 3 4

=-1 =-1/3 =1/3 =1

Figura 3.12 Elemento de cuatro nudos

N1 N2

+1
+1
 
1 2 3 4 1 2 3 4

Figura 3.13 Elemento de 4 nudos. Figura 3.14 Elemento de 4 nudos.

Interpolación de nudo esquina Interpolación de nudo interior
3.4. FORMULACIÓN ISOPARAMÉTRICA
En el apartado anterior se han definido las funciones de interpolación en la
coordenada local . Para calcular la matriz B es necesario hallar las derivadas de las
funciones de interpolación respecto de x, y por lo tanto se hace necesario establecer
algún tipo de relación entre la coordenada  y la x a fin de poder completar la
derivación.
En la formulación isoparamétrica, esta relación se introduce mediante las mismas
funciones de interpolación N usadas para interpolar la deformación, estableciendo la
relación:
x Ni xi N xe (3.28)

El vector xe contiene las coordenadas de todos los nudos del elemento.

Así, para el elemento de tres nudos, la ley de interpolación de las coordenadas es:
( 1) ( 1)
x x1 (1 )(1 ) x2 x3 (3.29)
2 2
La principal ventaja de esta formulación radica en que los tres nudos no tienen
porqué ser equidistantes, sino que el central (2) puede estar en cualquier posición
entre el 1 y el 3 (sin coincidir con ellos). En coordenadas locales, sin embargo, el nudo
central sigue estando en =0. Esta transformación de coordenadas corresponde a una
distorsión del elemento, que pasa a ser curvo, y cuyo mayor interés se manifiesta en
los elementos bidimensionales.

2

=-1 =0 =+1

x
x1 x2 x3

Figura 3.15 Interpolación de coordenadas

Es posible asimismo generar elementos subparamétricos, en los cuales las funciones

usadas para interpolar las coordenadas son de orden inferior a las usadas para
interpolar las deformaciones. Lo más habitual en este caso es usar funciones lineales
para interpolar las coordenadas, es decir apoyarse sólo en los dos nudos extremos
para definir la geometría, con lo que los nudos interiores deben estar centrados en el
elemento. El interés por los elementos subparamétricos se pone de manifiesto en el
problema bidimensional.
3.4.1 Interpolación de deformaciones
La ley de interpolación de deformaciones para un elemento general de n nudos es:

U1
e
u N 1 N 2 ... N n ... N (3.30)
Un

Empleando las funciones anteriores, esta ley será un polinomio de orden n-1.

3.4.2 Deformaciones unitarias. Matriz B

La expresión general de las deformaciones unitarias es:

U1
du e dN 1 dN n
N ... ... (3.31)
dx dx dx
Un

Esta expresión define la matriz B como:

dN 1 dN n
B ... (3.32)
dx dx
Para obtener las derivadas de las funciones de interpolación respecto de x, se aplica
la derivación en cadena:
dN i dN i d
(3.33)
dx d dx
La derivada de la función N respecto de la coordenada local  es conocida para cada
elemento. La derivada de una coordenada respecto de la otra se puede evaluar
empleando la ley de interpolación de coordenadas:
dx dN i
xi J (3.34)
d d
Donde se ha definido el jacobiano J de la transformación de coordenadas. Por lo
tanto queda:
dN i dN i 1
(3.35)
dx d J
La matriz B queda finalmente:
1 dN 1 dN n
B ... (3.36)
J d d
Obsérvese que esta matriz no puede calcularse si J=0, lo cual ocurre si la
transformación de coordenadas x/ no es correcta. Esto último puede ocurrir si
algunos de los nudos son coincidentes entre sí, o no están dispuestos de forma
ordenada en el elemento.

3.4.3 Matriz de rigidez

La expresión general de la matriz de rigidez es:

K BT D B dv (3.37)

En ella la matriz elástica es un escalar D=E, y el diferencial de volumen es:

dv Adx AJ d (3.38)

Sustituyendo la matriz B, se obtiene:

dN 1
1 d
1 dN 1 dN n E
K ... ... AJ d (3.39)
1
J d d J
dN n
d
Puede desarrollarse como:
dN 1 dN 1 dN 1 dN 2 dN 1 dN n
...
d d d d d d
dN 2 dN 1 dN 2 dN 2 dN 2 dN n
1
EA ...
K d d d d d d d (3.40)
1
J ... ... ... ...
dN n dN 1 dN n dN 2 dN n dN n
...
d d d d d d
Estudiando la naturaleza de los distintos términos del integrando se puede deducir
que si el jacobiano J es constante, el integrando es un polinomio. Sin embargo si J no
es constante, el integrando es un cociente de polinomios. En el primer caso la integral
puede efectuarse de forma exacta empleando métodos numéricos adecuados,
mientras que en el segundo la integración numérica siempre es aproximada.

3.4.4 Vector de fuerzas nodales equivalentes

Su expresión general, para el caso unidimensional es:

Pv NT q dv (3.41)
Donde q es la fuerza de volumen actuante sobre el elemento, que en este caso es un
escalar q, pues sólo tiene componente en x. Esta fuerza es en principio es variable,
pero con objeto de simplificar la práctica del método, es habitual restringir la posible
variación de la fuerza de volumen y limitarla sólo a aquellas que pueden ser
representadas por las propias funciones de interpolación. De esta forma la variación
de la fuerza de volumen se puede representar mediante una interpolación de sus
valores nodales, empleando las propias funciones de interpolación:
q N qev (3.42)

Siendo qev los valores nodales de la fuerza de volumen, que son constantes.

q3
q1 q2

1 q 2 3

Figura 3.16 Fuerza distribuida axial

Con esta aproximación, perfectamente válida a efectos prácticos, se obtiene la

siguiente expresión del vector de fuerzas nodales equivalentes:

Pv NT N dv qev M qev (3.43)

Considerando la estructura de la matriz N, la expresión de la matriz M es:

N 1N 1 N 1N 2 ... N 1N n
N 2N 1 N 2N 2 ... N 2N n
M NT N dv AJ d (3.44)
... ... ... ...
N n N 1 N n N 2 ... N n N n

El cálculo de la matriz M no presenta ningún problema, efectuándose en el sistema

local de coordenadas, en el que se conoce la matriz de interpolación N. Todos los
términos que la forman son polinomios, por lo que se puede evaluar de forma exacta
por métodos numéricos.

3.5. ELEMENTO ISOPARAMÉTRICO DE TRES NUDOS

La matriz de funciones de interpolación para este elemento es:
( 1) ( 1)
N (1 )(1 ) (3.45)
2 2
 U1 U2 U3

x1 x2 x3

Figura 3.17 Elemento isoparamétrico de tres nudos

El jacobiano de la transformación es:

dx dN i 2 1 2 1
J xi x1 2 x2 x3 (3.46)
d d 2 2
No es, en general, constante, y mide la distorsión en la transformación x/. Ocurren
varios casos:
 Si x1<x2<x3 el valor de J es siempre positivo.
 Si x1=x2, es decir que el elemento está colapsado, el valor de J es nulo en =-1/2.
 Si x2=x3, es decir que el elemento está colapsado, el valor de J es nulo en =+1/2.
 Si el nudo central está en el centro del elemento: x2=(x1+x3) , el valor de J es
constante:
x3 x1 L
J (3.47)
2 2
En este caso el elemento no está distorsionado
La matriz B es:
2 1 2 1
B 2 (3.48)
2 2
La matriz de rigidez vale:
2 1 2 2 1
2
1
EA 4 4
2 2
K 4 2 d (3.49)
1
J
2 1
simétrica
4
Para tratar elementos de área variable, se procede a interpolar ésta a través de sus
valores nodales, al igual que se ha hecho con las coordenadas:
A N Ae (3.50)

Siendo Ae los valores nodales del área. Por supuesto que esta aproximación limita la
variación del área a aquellos casos en los que dicha variación pueda representarse de
forma exacta mediante las funciones de interpolación.
La matriz M que proporciona las fuerzas nodales equivalentes vale:
2
( )2 ( 2
)(1 2
) ( 2
)( 2 )
4 2 4
1 2 2
2 2 ( )(1 )
M AJ (1 ) d (3.51)
1
2
2
( )2
simétrica
4
La distribución de tensiones es:
e E 2 1 2 1
E E B U1 2 U2 U3 (3.52)
J 2 2
Caso particular: elemento de tres nudos, de longitud L, con el nudo central centrado y
área constante A. Las expresiones anteriores se simplifican pues J L / 2 .
Matriz de rigidez:

7 8 1
EA
K 8 16 8 (3.53)
3L
1 8 7

Matriz de fuerzas nodales equivalentes:

4 2 1
AL
M 2 16 2 (3.54)
30
1 2 4

Distribución de tensiones: es lineal en el elemento, por ser J constante.

E
(2 1)U 1 4 U2 (2 1)U 3 (3.55)
L
 4
Elasticidad bidimensional

4.1. INTRODUCCIÓN
Los problemas de elasticidad bidimensional son muy frecuentes en Ingeniería, y son
asimismo los primeros en los que se aplicó el MEF. En este caso el medio continuo
que se analiza es plano, y se considera situado en el plano XY. Se denomina t al
espesor del dominio en su dirección transversal, el cual se considera despreciable
frente a las dimensiones del dominio en el plano XY.
La posición de un punto está definida por dos coordenadas (x,y), y su deformación
tiene dos componentes u(x,y), v(x,y) en las direcciones x,y respectivamente. El
campo de deformaciones es por lo tanto un vector:
u(x , y )
u (4.1)
v(x , y )

Dentro de la elasticidad en dos dimensiones existen dos problemas diferentes:

 Tensión plana: cuando la tensión z en sentido perpendicular al plano XY es cero,
ya que el sólido puede dilatarse libremente en el sentido de su espesor. Por lo
tanto existe una deformación unitaria z no nula en dicha dirección.
 Deformación plana: cuando en el sentido del espesor del sólido no hay posibilidad
de deformación, es decir z=0 por lo que se genera una tensión en dicha dirección
z no nula.
En ambos casos la tensión y la deformación en la dirección z no contribuyen a la
energía elástica del sistema.
Las ecuaciones diferenciales que rigen el problema son de orden m=2 en las
deformaciones, pues contienen la derivada primera de las tensiones, que a su vez son
derivadas de las deformaciones. En la expresión del potencial total del sistema (ver
apartado 2.7) aparecen las deformaciones unitarias , que son las derivadas primeras
de las deformaciones u, luego el orden de derivación de las incógnitas primarias u en
el potencial es n=1. Se requieren por lo tanto funciones de interpolación con
continuidad C0 para asegurar la convergencia del método.

4.2. FUNCIONES DE INTERPOLACIÓN

El campo de deformaciones en el interior del elemento se aproxima mediante la
expresión habitual:
u Ni Ui v NV
i i
(4.2)

En forma matricial es:

e
u N (4.3)
e
Las deformaciones nodales del elemento se agrupan en un vector columna:
T
e
U 1 V1 U 2 V2 ... U n Vn (4.4)

siendo n el número de nudos del elemento. La matriz de funciones de interpolación

N tiene 2 filas y tantas columnas como grados de libertad haya entre todos los nudos
del elemento. La estructura de esta matriz siempre es la misma:
N1 0 N2 0 ... 0 Nn 0
N (4.5)
0 N1 0 N2 0 ... 0 Nn

V1
U1

v
V2 u
U2
V3

U3

Figura 4.1 Interpolación de deformaciones

4.3. DEFORMACIONES UNITARIAS

Las deformaciones unitarias en un punto del elemento finito son:
 u
x
x
v
y (4.6)
y
xy
u v
y x
Se pueden poner en la forma:

0
x
x u
0 u (4.7)
y v
y

xy

y x
donde se identifica al operador matricial  que pasa de las deformaciones u a las
deformaciones unitarias. Sustituyendo las deformaciones u en función de las
deformaciones nodales, a través de las funciones de interpolación, se obtiene:
e e
u N B (4.8)

Se identifica de esta forma la matriz B:

N1 N2 Nn
0 0 ... 0
x x x
N1 N2 Nn
B N 0 0 ... 0 (4.9)
y y y
N1 N1 N2 N2 Nn Nn
...
y x y x y x
Esta matriz relaciona las deformaciones de los nudos con las deformaciones unitarias
en un punto cualquiera del elemento. Por lo tanto B representa el campo de
deformaciones unitarias que se supone existe en el interior del elemento finito, como
consecuencia de la hipótesis de interpolación de deformaciones efectuada.
Esta matriz se puede poner en la forma:

B B1 B2 ... Bn (4.10)

Siendo cada una de las submatrices Bi

 Ni
0
x
Ni
Bi 0 (4.11)
y
Ni Ni
y x
Nótese que debido a la estructura de B, las deformaciones unitarias en el interior del
elemento se pueden poner en función de las deformaciones nodales en la forma:
Ni Ni Ni Ni
x U y V xy U V (4.12)
x i y i y i x i

4.4. ESTADO DE TENSIONES. ECUACIÓN CONSTITUTIVA

El estado de tensiones en dos dimensiones es:

y (4.13)
xy

La ecuación constitutiva, en ausencia de temperaturas, es:

D (4.14)

Para un material elástico lineal e isótropo la matriz elástica D es constante. Su

expresión es diferente para los dos problemas de elasticidad plana.
 Tensión plana. En este caso la tensión en dirección transversal al material (z) es
nula, pero éste es libre de dilatarse en dicha dirección z: z
0 z
0

Se parte de la ecuación constitutiva en el estado tridimensional:

2 0 0 0
x x

y
2 0 0 0 y

z 2 0 0 0 z

xy 0 0 0 0 0 xy
(4.15)
yz 0 0 0 0 0 yz

zx
0 0 0 0 0 zx

E E
(1 )(1 2 ) 2(1 )
Imponiendo en la tercera ecuación la condición z 0 se obtiene:
x y ( 2 ) z 0 (4.16)

De donde se calcula el valor de la deformación unitaria transversal al material:

z ( x y ) (4.17)
2
Sustituyendo en la expresión inicial del estado tridimensional (y considerando
además que yz=0, zx=0), se obtiene la matriz elástica del estado plano:

1 0
E
DTP 1 0 (4.18)
1 v2
1
0 0
2

y
xy
x x

Y z
z=0
y
Z X

Figura 4.2 Estado de tensión plana

 Deformación plana. En este caso la deformación unitaria transversal al material

(z) es nula, pues éste es incapaz de dilatarse en dirección z. En consecuencia debe
existir tensión en dicha dirección, es decir: z
0 z
0

Para obtener la ecuación constitutiva es suficiente con hacer cero las deformaciones
unitarias correspondientes en la ecuación tridimensional: basta por lo tanto con
extraer las filas y columnas correspondientes al estado plano. De esta forma se
obtiene la siguiente matriz elástica:

1 0
1
E (1 )
DDP 1 0 (4.19)
(1 )(1 2 ) 1
1 2
0 0
2 2
Nótese que aparece una tensión en la dirección z, cuyo valor se deduce simplemente
de la ecuación en la dirección z:
E
z ( x y ) (4.20)
(1 )(1 2 )

4.5. DEFORMACIONES UNITARIAS INICIALES. TEMPERATURAS

La presencia de deformaciones unitarias iniciales introduce un nuevo término en la
ecuación constitutiva:
D( 0
) (4.21)

siendo 0 el vector de las deformaciones unitarias iniciales existentes en el material

en el punto considerado, que deben ser conocidas.
Si las deformaciones unitarias iniciales 0 están producidas por un incremento de
temperatura T, su valor es, para el estado de tensión plana:
T
0Tx

0 0Ty
T (4.22)
0Txy 0

Las deformaciones unitarias iniciales son iguales en ambas direcciones x e y, y

además no se genera deformación de cortante.
Para el caso de deformación plana, a la dilatación T en las dos direcciones x,y se
suma la imposibilidad de dilatarse según z, lo que origina una tensión en dicho
sentido, la cual genera a su vez una deformación unitaria según los ejes x,y, por
efecto de Poisson, y que se suma a las iniciales:
T
0Tx

0 0Ty
(1 ) T (4.23)
0Txy 0

4.6. ELEMENTO TRIANGULAR

Este elemento tiene seis desplazamientos en los nudos, que forman un vector:
T
e
U 1 V1 U 2 V2 U 3 V3 (4.24)

Los desplazamientos de un punto cualquiera dentro del elemento se pueden

representar en función de estos seis valores nodales, mediante una expresión
polinómica. Dado que hay seis deformaciones nodales, el polinomio sólo podrá tener
seis términos:
u 1
x
2 3
y v 4
x
5
y
6
(4.25)

Estas expresiones pueden ponerse en forma matricial:

u 1 x y 0 0 0 1

v ... (4.26)
0 0 0 1 x y
6

u R (4.27)

Los seis parámetros i se pueden calcular aplicando la expresión (4.26) a los tres
nudos del elemento, y agrupando las seis ecuaciones obtenidas (dos en cada nudo):
U1 1 x1 y1 0 0 0
1
V1 0 0 0 1 x1 y1 2

U2 1 x2 y2 0 0 0 3
(4.28)
V2 0 0 0 1 x2 y2 4

U3 1 x3 y3 0 0 0 5

V3 0 0 0 1 x 3 y3 6

es decir:
e
C (4.29)

Despejando  y sustituyendo en la expresión de u se obtiene:

1 e
u RC (4.30)

Esta expresión define las funciones de interpolación como:

1
N RC (4.31)

Efectuando el producto de matrices anterior se obtiene la expresión:

N1 0 N2 0 N3 0
N (4.32)
0 N1 0 N2 0 N3

Las tres funciones de interpolación correspondientes a los tres nudos son:

a1 b1x c1y a2 b2x c2y a3 b3x c3y
N1 N2 N3 (4.33)
2A 2A 2A

Las distintas constantes dependen de la geometría del elemento:

 a1 x 2y3 x 3y2 b1 y2 y3 c1 x3 x2

a2 x 3y1 x1y3 b2 y3 y1 c2 x1 x3 (4.34)

a3 x1y2 x 2y1 b3 y1 y2 c3 x2 x1

En la ecuación anterior, A es el área del elemento, cuyo valor se obtiene mediante el

determinante:

1 x1 y1
1
A Det 1 x 2 y2 (4.35)
2
1 x 3 y3

Se observa que si el elemento tiene área nula (dos nudos coincidentes) eso se
manifiesta en A=0 y no se pueden calcular las Ni. Estas funciones son planos de valor
1 en el nudo i y 0 en los otros dos nudos.
Al ser las tres funciones de interpolación planos (Figura 4.3), el campo de
desplazamientos en el interior del elemento es también un plano que pasa por los
tres valores nodales del desplazamiento. En consecuencia, si se emplea este
elemento, el campo de deformaciones en una dirección cualquiera u o v se aproxima
mediante una superficie poliédrica de facetas triangulares.

N1 N2 N3
1
2 1 2 1
2

3
3 3

Figura 4.3 Funciones de interpolación del elemento triangular

 El estado de deformación unitaria viene definido por la matriz B, que en este caso
vale:
N1 N2 N3
0 0 0
x x x
N1 N2 N3
B N 0 0 0 (4.36)
y y y
N1 N1 N2 N2 N3 N3
y x y x y x
Efectuando las derivadas, la expresión de esta matriz es:
 b1 0 b2 0 b3 0
1
B 0 c1 0 c2 0 c3 (4.37)
2A
c1 b1 c2 b2 c3 b3

Esta matriz se puede poner en la forma:

B B1 B2 B3 (4.38)

Se observa que la matriz B es constante, y no depende de x,y. Por lo tanto las

deformaciones unitarias  son constantes en todo el elemento, y también lo serán las
tensiones, que son proporcionales a ellas.
 La matriz de rigidez del elemento se obtiene a partir de su expresión general.
Dada la estructura de B, se puede poner en la forma:

BT1
K BT D B dv BT2 D B1 B2 B3 t dxdy (4.39)
v
BT3

siendo t el espesor del elemento. La matriz K se puede dividir en 3x3 submatrices,

que relacionan a los tres nudos entre sí:

K11 K12 K13

K K21 K22 K23 (4.40)
K31 K32 K33

Cada una de ellas tiene la expresión:

Kij BTi D Bj t dxdy (4.41)

Dado que Bi y D son constantes, su valor resulta ser (suponiendo espesor t uniforme):

ij t bibj D11 cic j D33 bic j D12 cibj D33

K (4.42)
4A cibj D12 bic j D33 cic j D11 bibj D33

donde Dij son los coeficientes de la matriz elástica D.

 El estado de tensiones en un punto cualquiera del elemento es:
e
D DB (4.43)
 U1
1 0 b1 0 b2 0 b3 0
x
E 1 V1
1 0 0 c1 0 c2 0 c3 (4.44)
y
(1 2
) 2A ...
xy 1 c1 b1 c2 b2 c3 b3
0 0 V3
2
Se observa que el estado de tensiones no depende de las coordenadas x,y sino que
es uniforme en todo el elemento. Por lo tanto, si se emplea este elemento el campo
de tensiones en el material se aproxima mediante una serie de valores constantes,
formando una superficie escalonada. En consecuencia este elemento sólo debe
utilizarse con mallados muy finos si se prevé un campo de tensiones variable.
 Las fuerzas nodales equivalentes a las deformaciones iniciales térmicas son:

BT1
PT BT D 0T
dv BT2 D 0T
t dx dy (4.45)
v
BT3

Suponiendo un incremento de temperatura T uniforme en todo el elemento, el

integrando de la expresión anterior resulta ser constante con lo que se obtiene, para
el caso de tensión plana:
PT 1x b1
PT 1y c1
PT 2x EA Tt b2
PT (4.46)
PT 2y 1 c2
PT 3x b3
PT 3y c3

4.7. ELEMENTO RECTANGULAR

El desarrollo de este elemento sigue un camino paralelo al del triángulo. Para mayor
sencillez en la formulación se utiliza en este elemento el sistema de ejes local
mostrado en la figura 4.4.
 El vector columna de desplazamientos nodales es:
T
e
U 1 V1 U 2 V2 U 3 V3 U 4 V4 (4.47)
 V2 V1 YL
U2 U1 2 1
v
u XL

2c
V3 V4
U3 U4
3 4
2b

Figura 4.4 Elemento rectangular

Con los ocho grados de libertad, el campo de desplazamientos en el interior del

elemento se puede representar mediante dos polinomios bilineales, de cuatro
términos:

u 1 x y xy 0 0 0 0 1

v ... (4.48)
0 0 o 0 1 x y xy
8

u R (4.49)

El cálculo de  se efectúa nuevamente particularizando la ecuación anterior para los

cuatro nudos, a base de sustituir en x,y las coordenadas de dichos nudos, que con el
sistema de ejes elegido son muy sencillas. Queda una ecuación del tipo:
e
C (4.50)

La matriz C es de tamaño 8x8 y sus términos son todos conocidos en función del
tamaño del elemento. Por ejemplo sus dos primeras filas son:
1 b c bc 0 0 0 0
0 0 0 0 1 b c bc
C (4.51)
1 b c bc 0 0 0 0
.. .. .. .. .. .. .. ..
1 e
Despejando C y sustituyendo en la ecuación (4.49) se obtiene:
1 e
u R RC (4.52)

Se identifica la matriz de funciones de interpolación N R C 1 . Operando se llega a:

N1 0 N2 0 N3 0 N4 0
N (4.53)
0 N1 0 N2 0 N3 0 N4

Las distintas funciones de interpolación son del tipo bilineal:

 N1 (b x )(c y )/ A
N2 (b x )(c y )/ A
(4.54)
N3 (b x )(c y )/ A
N4 (b x )(c y )/ A

La figura 4.5 muestra el aspecto de dos de estas funciones.

N1 N2
2 2
1 1

3 4 3 4

Figura 4.5 Funciones de interpolación del elemento rectangular

 La matriz B para este elemento es de 3x8, y está formada por cuatro matrices Bi
de 3x2, una para cada nudo:

B B1 B2 B3 B4 (4.55)

(c y) 0 (c y) 0 (c y) 0 (c y) 0
1
B 0 (b x) 0 (b x) 0 (b x) 0 (b x)
A
(b x ) (c y) (b x) (c y) (b x) (c y) (b x) (c y)

Se observa que esta matriz tiene términos lineales, por lo que ésta es la variación
permitida para el campo de deformaciones unitarias en el interior del elemento. De
la misma forma, las tensiones también pueden tener una variación lineal en el
elemento. Este elemento es por lo tanto bastante más preciso que el triangular, que
sólo permitía valores constantes de la tensión y la deformación unitaria.
 La matriz de rigidez de este elemento es de tamaño 8x8, y se calcula utilizando la
expresión habitual:

BT1
K BT D B dv ... D B1 ... B4 t dxdy (4.56)
BT4

La matriz K se puede dividir en 4x4 submatrices, que relacionan a los cuatro nudos
entre sí. Cada una de ellas vale:

Kij BTi D Bj t dxdy (4.57)

Esta matriz se obtiene referida al sistema de ejes local del elemento. Para pasarla al
sistema general se aplica una transformación de coordenadas plana.

4.8. FUNCIONES DE INTERPOLACIÓN DE TIPO C0

En los apartados anteriores se han desarrollado elementos finitos sencillos, a base de
definir un polinomio interpolante de las deformaciones y a continuación determinar
los coeficientes de dicho polinomio a base de ajustarlo a los valores nodales de las
deformaciones. Una vez ajustado este polinomio y determinadas las funciones de
interpolación Ni, el proceso de cálculo de las propiedades del elemento es siempre el
mismo, con independencia del tipo de elemento. La definición de las funciones de
interpolación es por lo tanto el paso fundamental en el análisis por el MEF, y de él
depende en gran manera la precisión de los resultados obtenidos.
La matriz N tiene tantas filas como desplazamientos se consideren para un punto del
continuo (dos en el caso plano), y tantas columnas como grados de libertad haya
entre todos los nudos del elemento. La estructura de esta matriz siempre es la
misma, y para el caso plano es:
N1 0 N2 0 ... 0 Nn 0
N (4.58)
0 N1 0 N2 0 ... 0 Nn

siendo n el número de nudos del elemento. Resulta por lo tanto de gran interés
definir las funciones de interpolación de los nudos Ni, para cada tipo de elemento.
En los apartados siguientes se presentan funciones de interpolación con
compatibilidad de tipo C0, es decir que garantizan la continuidad de la propia función
en los bordes entre elementos.
Las funciones de interpolación siempre son del tipo polinómico. El número de
términos de éste polinomio viene determinado por el número de grados de libertad
del elemento, que define el número de parámetros independientes que pueden
utilizarse para definir el polinomio. En general se trata de utilizar polinomios
completos del mayor grado posible. El número de términos que aparecen en un
polinomio de grado dado y dos variables x,y, se puede deducir del triángulo de
Pascal. Por ejemplo: un polinomio completo de grado 1 requiere tres términos, el de
grado 2 requiere seis términos, el de grado 3, diez términos, etc.

4.8.1 Elementos rectangulares. Formulación Serendipity

En esta formulación las funciones de interpolación se definen de la siguiente manera:
 Se utilizan únicamente nudos situados en los contornos del elemento.
 Se utiliza un sistema de coordenadas local al elemento , con origen en su
centro y normalizadas de tal forma que valen +1 y -1 en sus extremos. En este
 sistema de coordenadas local, el elemento queda representado como un
cuadrado de lado 2.
 El número de nudos existente en cada lado del elemento define el grado del
polinomio interpolado en dicho lado.
 La función de interpolación completa es el producto de las funciones en cada
dirección.


(+1,+1)

(-1,-1)

Figura 4.6 Coordenadas locales normalizadas

4.8.1.1 Elemento de cuatro nudos

Las funciones de este elemento son:
1
Ni (1 i
)(1 i
) i 1, 4 (4.59)
4
Siendo i,i, las coordenadas del nudo i. Esta funciones son bilineales, y son las
mismas que se dedujeron para el elemento rectangular, aunque ahora están
expresadas en las coordenadas locales . Para desarrollos posteriores, a las
funciones de este elemento se les denomina funciones básicas, y se designan como
N̂ .


2 1

3 4

Figura 4.7 Elemento de cuatro nudos

4.8.1.2 Elemento de ocho nudos

En cada lado las funciones de interpolación varían cuadráticamente, ya que hay tres
nudos en cada uno. Adoptando la numeración de la figura 4.8, las funciones para los
nudos centrales de los distintos lados son:
 1 2
Nk (1 )(1 k
) k 2 6 k
0
2
1 2
Nk (1 )(1 k
) k 4 8 k
0 (4.60)
2


2 1
3

8

4

5 6 7

Figura 4.8 Elemento de ocho nudos

6
3

2
7

8
1

Figura 4.9 Elemento de 8 nudos. Función de interpolación de un nudo central

Al introducirse un nudo en el centro de un lado, las funciones de los nudos de las

esquinas se calculan modificando las básicas Nˆi de la forma siguiente (figura 4.10):

Ni Nˆi 1
2 Nk1 1
2 Nk 2 (4.61)

k2 i

k1

Figura 4.10 Corrección de nudos esquina

Es decir que a cada función básica de una esquina se le resta 1/2 de la función de los
dos nudos intermedios de los lados adyacentes. Efectuando dicha operación para los
cuatro nudos de las esquinas se obtienen sus funciones de interpolación:
 1
Ni (1 i
)(1 i
)( i i
1) i 1 3 5 7 (4.62)
4

5
4

7 2

8
1

Figura 4.11 Elemento de 8 nudos. Función de interpolación de nudo esquina

4.8.1.3 Elemento de doce nudos

En cada lado del elemento la variación de las funciones de interpolación es cúbica, al
disponerse de cuatro nudos para definirla.

3 2
4 1

5 12 

6 11

7 8 9 10

Figura 4.12 Elemento de doce nudos

Los dos nudos intermedios de un lado cualquiera corresponden a valores +1/3 y -1/3
de la coordenada. Las funciones de dichos nudos intermedios se pueden obtener de
forma intuitiva, definiendo una cúbica que adopte los valores deseados en los nudos:
Nk 9 (1 )(1 9 )(1 2
) i 5 6 11 12 1
32 k k k

Nk 9 (1 )(1 9 )(1 2
) i 2 3 8 9 1 (4.63)
32 k k k

Para los nudos de esquina se parte de las funciones básicas del elemento de cuatro
nudos Nˆi . Al introducir los nudos intermedios las funciones básicas se ven afectadas
de la forma siguiente:
Ni Nˆi 2
3 Nk1 1
3 Nk 2 (4.64)
donde i es el nudo esquina, y k1 y k2 son los dos nudos intermedios adyacentes,
siendo k1 el más próximo a i y k2 el más lejano (figura 4.13).
Para el otro nudo esquina del lado, que llamaremos j, es:
Nj Nˆ j 2
3 Nk 2 1
3 Nk1 (4.65)

Esta corrección debe hacerse para todos los lados y todas las esquinas del elemento.
Se obtiene la expresión siguiente, que es válida para los cuatro nudos esquina:
2
Ni 1
32 (1 i
)(1 i
)( 10 9 9 2) i 1 4 7 10 (4.66)

j k2 k1 i

Figura 4.13 Corrección de nudos esquina

4.8.1.4 Elementos con número de nudos variable

Estos elementos se usan para servir de transición entre zonas discretizadas con
distintos tipos de elementos, como se muestra en la figura 4.14.
La obtención de sus funciones de interpolación puede hacerse de forma directa, pero
dada la gran cantidad de casos que pueden presentarse, lo mejor es seguir un camino
sistemático:
 Para los nudos intermedios de los lados, se calculan sus Nk según las expresiones
antes indicadas, en función de cuántos nudos intermedios haya en ese lado.
 Para los nudos de las esquinas se corrigen las funciones básicas del elemento de
cuatro nudos Nˆi , en función de los nudos que aparezcan en los lados adyacentes
a cada esquina, como se ha indicado para los elementos de tres y cuatro nudos
por lado.

8 nudos Transición 4 nudos

Figura 4.14 Elementos con número de nudos variable

4.8.1.5 Grado de los polinomios generados
Los términos de los polinomios utilizados por los distintos elementos en formulación
Serendipity se pueden observar en el triángulo de Pascal:

1
4 nudos
x y
8 nudos
x2 xy y2
12 nudos
x3 x2y xy2 y3
x4 x3y x2y2 xy3 y4
x5 x4y x3y2 x2y3 xy4 y5
x6 x5y x4y2 x3y3 x2y4 xy5 y6

Se observa que disminuye el número de términos de grado elevado, con respecto a la

interpolación de Lagrange. Por ejemplo en el elemento cuadrático se disminuye el
número de términos en exceso de tres a dos, y en el cúbico se disminuye de seis a
dos.
Sin embargo utilizando únicamente nudos situados en los contornos, como máximo
se puede generar un polinomio completo hasta grado 3. Si se quiere pasar a grado
cuatro o superior, hay que añadir nudos intermedios. Esto ocurre en el elemento de
cinco nudos por lado, que tiene 17 nudos (16 en los lados y uno interior).

4.8.2 Elementos rectangulares - Interpolación de Lagrange

Se utilizan polinomios de Lagrange como funciones de interpolación, que tienen la
propiedad de valer 1 en un nudo y 0 en los demás nudos. Se utilizan las coordenadas
locales al elemento ,, normalizadas de tal forma que valen +1 y -1 en los extremos
del elemento.
En una dimensión, el polinomio de Lagrange de grado n que pasa por el nudo i, se
ajusta a n+1 puntos, y vale 1 en el nudo i, y cero en los demás (Figura 4.15):
( j
)
lin ( ) (4.67)
j 0,n ; j i ( i j
)

4
l1
Figura 4.15 Interpolación de Lagrange

En dos dimensiones, la función de interpolación del nudo i se obtiene como producto

de dos polinomios de Lagrange, uno en cada dirección. El orden de cada polinomio
depende del número de nudos en cada dirección:
 ( j
) ( )
Ni ( , ) lin ( )lim ( ) k
(4.68)
j 0,n ; j i ( i j
)k 0,n ;k i ( i k
)

siendo n y m el número de divisiones en cada dirección.

Figura 4.16 Elementos con interpolación de Lagrange en dos dimensiones

Estas funciones Ni son muy fáciles de generar, pero tienen dos problemas: el primero
es el gran número de nudos que introducen en el elemento, aumentando el número
de ecuaciones a resolver. El segundo es que no utilizan polinomios completos, sino
que consideran muchos términos de grado elevado, que son peores para ajustar la
solución, pues introducen ondulaciones en ella (términos parásitos). Sin embargo,
por esta característica, son más adecuados para tratar problemas de transmisión de
ondas que requieren de términos de orden alto en la solución.

1
4 nudos
x y
9 nudos
x2 xy y2
16 nudos
x3 x2y xy2 y3
x4 x3y x2y2 xy3 y4
x5 x4y x3y2 x2y3 xy4 y5
x6 x5y x4y2 x3y3 x2y4 xy5 y6

4.8.3 Elementos triangulares

En los elementos triangulares se utilizan unas coordenadas locales muy particulares,
que son las coordenadas de área (figura 4.17). Estas coordenadas se definen como:
Area(P 23)
L1 (4.69)
Area(123)
y análogamente las otras 2 coordenadas L2 y L3.
La expresión general es:
Area(PJK )
LI (4.70)
Area(IJK )
Evidentemente se cumple que L1 + L2 + L3 = 1
Las líneas de Li constante son rectas paralelas al lado JK (figura 4.18).

1 L1=1 1

0.6
L3 L2 0.4
0.2
L1
2 3 L1=0 3
2

Figura 4.17 Coordenadas de área Figura 4.18 Valores constantes de la

coordenada de área L1

4.8.3.1 Elemento de tres nudos

Las funciones de interpolación son: N i Li i 1, 3
La función del nudo i es un plano que pasa por el lado opuesto al vértice i, y que vale
1 en dicho nudo. Al igual que en los rectángulos estas funciones se llaman básicas
Nˆi .
1

2 3

Figura 4.19 Elemento de tres nudos

4.8.3.2 Elemento de seis nudos

En sus lados las funciones de interpolación varían de forma cuadrática. Cada nudo
intermedio tiene una función que vale:
N2 4L1L2 N4 4L2L3 N6 4L3L1 (4.71)
 N2
1

2 6
2

3 5 5
4 3 4

Figura 4.20 Elemento de seis nudos y función de interpolación del nudo 2

Las funciones de los nudos de esquina se calculan a partir de las funciones básicas,
restando a cada una de ellas 1/2 de las funciones correspondientes a los nudos
intermedios que hay en los lados adyacentes (figura 4.21):
Ni Nˆi 1
2 Nk1 1
2 Nk 2 (4.72)

k1 k2

Figura 4.21 Corrección de nudo esquina

Se obtienen las siguientes funciones de interpolación:

N1 (2L1 1)L1 N3 (2L2 1)L2 N5 (2L3 1)L3 (4.73)

4.8.3.3 Elemento de diez nudos

En este elemento se debe incluir un nudo en su centro con el fin de cubrir todo un
polinomio completo de grado tres, que requiere diez términos.

2 9

3 10
8

4 7
5 6

Figura 4.22 Elemento triangular de diez nudos

Las funciones de los nudos intermedios son:
N2 9
2 L1L2 (3L1 1) N 3 9
2 L1L2 (3L2 1)
N5 9
2 L2L3 (3L2 1) N 6 9
2 L2L3 (3L3 1) (4.74)
N8 9
2 L3L1(3L3 1) N 9 9
2 L3L1(3L1 1)

Para el nudo central la función de interpolación es la función de la burbuja:

N10 27L1L2L3 , que tiene valor nulo en todo el contorno del triángulo.
Las funciones de las esquinas se obtienen de manera análoga al caso de los
rectángulos. Cada lado con dos nudos afecta a los dos nudos de esquina del lado de
la forma siguiente (figura 4.23):
Ni Nˆi 2
3 Nk1 1
3 Nk 2 (4.75)
donde i es el nudo esquina, y k1 y k2 son los dos nudos intermedios adyacentes,
siendo k1 el más próximo a i y k2 el más lejano. Para el otro nudo esquina del lado,
que llamaremos j la corrección es:
Nj Nˆ j 2
3 Nk 2 1
3 Nk1 (4.76)

j k2 k1 i

Figura 4.23 Corrección de nudos esquina

4.9. FORMULACIÓN ISOPARAMÉTRICA

Se han estudiado hasta el momento diversos tipos de elementos que se caracterizan
por un número creciente de nudos, lo que les permite interpolar el campo de
desplazamientos de forma sofisticada, y alcanzar una gran precisión en el análisis. Sin
embargo para todos ellos se han utilizado formas sencillas como triángulos o
rectángulos, mediante los cuales resulta difícil el adaptarse a la forma de la
estructura a estudiar, que normalmente tiene forma curva. En este apartado se trata
de estudiar elementos finitos con forma más compleja: así, por ejemplo un elemento
que es rectangular en sus coordenadas locales puede ser transformado en uno que
es curvilíneo en coordenadas cartesianas.
Para ello es necesario definir alguna relación entre las coordenadas de un punto en el
sistema local y las mismas en el sistema general:
x f( , ) y g( , ) (4.77)
Esta transformación de coordenadas define la forma del elemento en el sistema
cartesiano.

4.9.1 Interpolación de coordenadas

Uno de los métodos más convenientes para definir la relación entre coordenadas
locales y generales anterior consiste en usar funciones de interpolación, de la misma
forma que se hizo para el campo de deformaciones dentro del elemento. Por ejemplo
se puede poner:
x N i'x i y N i'yi (4.78)
1,n 1,n

donde N' son algunas funciones de interpolación definidas en las coordenadas

locales ,. Estas N' deben de valer 1 en el nudo i y cero en los demás. De inmediato
se ve la posibilidad de utilizar para estas N' a las funciones de interpolación
empleadas para interpolar los desplazamientos, que satisfacen dicha condición.

 

 N
2

2 X

Figura 4.24 Interpolación de coordenadas

En principio no hay por qué utilizar como funciones de interpolación de coordenadas

a las mismas funciones de interpolación de desplazamientos, sino que pueden ser
otras. A los elementos que utilizan las mismas funciones para interpolación de
coordenadas que las usadas para la interpolación de desplazamientos, se les llama
isoparamétricos. Es decir que estos elementos usan las coordenadas de todos sus
nudos para definir el cambio de coordenadas, o lo que es lo mismo, la forma del
elemento.
Por lo tanto, en la formulación isoparamétrica se puede poner:
x1
y
x N 1 0 N 2 0 ... 1
x2 (4.79)
y 0 N 1 0 N 2 ...
y2
...
 x N xe (4.80)

El vector xe agrupa a las coordenadas (x,y) de todos los nudos del elemento.
Con esta definición un elemento isoparamétrico tiene una forma geométrica real que
está definida por el tipo de funciones de interpolación que emplea. La forma real de
cada lado está definida por el número de nudos que hay en ese lado: así los lados con
dos nudos son rectos (funciones lineales), los lados de tres nudos pueden ser
parábolas (las funciones son cuadráticas) y los lados de cuatro nudos pueden ser
curvas cúbicas. En el sistema local  del elemento, éste siempre es un cuadrado de
lado 2.

Figura 4.25 Elementos con lados curvos

Elementos sub y super paramétricos

Otro tipo de elementos que pueden definirse son los subparamétricos. En ellos se
utilizan como funciones N' para interpolar las coordenadas, una funciones de orden
inferior a las de interpolación de desplazamientos. Es decir que no se utilizan las
coordenadas de todos los nudos para definir la forma del elemento, sino únicamente
las de algunos de ellos (normalmente las esquinas).
Así, un elemento de tres nudos por lado (Ni parabólicas) puede utilizar sólo los nudos
esquina para definir las coordenadas, usando unas N' lineales: en este caso sus lados
son forzosamente rectos, y el nudo central está siempre en el centro del lado (figura
4.26). Análogamente ocurre para el elemento de cuatro nudos por lado, en que se
emplean los 4 nudos de esquina para definir sus coordenadas.
Estos elementos son prácticos cuando la geometría es sencilla, y el elemento tiene
lados rectos, pero el campo de deformaciones varía mucho y debe aproximarse con
funciones de orden superior. Computacionalmente no tienen una gran ventaja, pero
permiten ahorrar gran cantidad de datos en la definición de las coordenadas de los
nudos, pues basta definir las coordenadas de las esquinas.
Los elementos superparamétricos son aquellos que utilizan más puntos para definir
su forma geométrica, que para definir el campo de deformaciones. Es decir que las N'
son de orden superior a las N. Así por ejemplo se pueden usar los elementos de la
figura 4.27.
 Deformaciones
Coordenadas

Figura 4.26 Elemento subparamétrico

Deformaciones
Coordenadas

Figura 4.27 Elemento superparamétrico

El primero permite una forma parabólica del contorno del elemento y sólo una
variación lineal en los desplazamientos. Estos elementos no son usados en la
práctica, pues requieren dar gran cantidad de coordenadas de nudos y sin embargo
no aprovechan dichos nudos para interpolar los desplazamientos.

4.9.2 Matriz de rigidez

La expresión obtenida para la matriz de rigidez de un elemento plano cualquiera es:

K BT D B dv (4.81)
v

En ella aparece la matriz B, que se obtiene derivando la matriz de funciones de

interpolación N respecto a las coordenadas x e y. Sin embargo, la matriz N está
definida en las coordenadas locales normalizadas del elemento , por lo que es
necesario transformar las derivadas entre unas coordenadas y otras. Una vez
obtenida la matriz B hay que efectuar la integral anterior al dominio de todo el
elemento. Ésta es muy simple en coordenadas locales pero es complicada en general
en coordenadas generales, dado que el elemento puede ser curvo.
Para resolver este problema lo que se hace es evaluar la integral anterior en
coordenadas locales , pasando todas las expresiones necesarias a dichas
coordenadas. En la expresión de la matriz B se observa que es necesario disponer de
todas las derivadas de las funciones de interpolación respecto a las coordenadas
generales x,y.
Siguiendo las reglas de la derivación se puede poner la siguiente relación matricial
entre las derivadas en ambos sistemas de coordenadas:
Ni x y Ni Ni
x x
J (4.82)
Ni x y Ni Ni
y y

La matriz J es la Jacobiana de la transformación de coordenadas (x,y) a .

Despejando se obtiene:
Ni Ni
x
J 1
(4.83)
Ni Ni
y

El vector de la derecha es conocido sin más que derivar las Ni respecto a , y
conociendo J se pueden obtener de la expresión anterior todas las derivadas que
forman la matriz Bi.
 El cálculo de J se hace apoyándose en la interpolación de coordenadas:
Ni Ni
xi yi
J (4.84)
Ni Ni
xi yi

Esta expresión puede ser evaluada fácilmente ya que las funciones N son conocidas
en función de  y xi,yi son las coordenadas de los nudos que definen la forma del
elemento.
 El dominio de integración expresado en coordenadas locales es:
dv t dxdy tJ d d (4.85)

donde aparece el determinante J de la matriz Jacobiana, que se calcula a partir de la

expresión de ésta, y es en general una función de .
La presencia del espesor del elemento t en la integral anterior permite tratar con
gran sencillez elementos de espesor variable. Para ello la solución práctica más
habitual consiste en suponer que la ley de variación del espesor puede aproximarse
por interpolación de los valores del espesor en los nudos ti, mediante la expresión
habitual:
t N iti (4.86)

 Considerando la estructura de la matriz B, la matriz de rigidez K es:

BT1
K BT D B dv ... D B1 ... Bn t dxdy (4.87)
BTn

La matriz K se puede dividir en n x n submatrices, que relacionan a los n nudos entre

sí:

K11 ... K1n

K ... ... ... (4.88)
Kn 1 ... Knn

Cada una de las submatrices tiene la expresión:

1
ij
K BTi D B j t J d d (4.89)
1

Sustituyendo los valores de las distintas matrices se obtiene:

Ni Nj Ni Nj Ni Nj Ni Nj
1 D11 D33 D12 D33
x x y y x y y x
Kij tJ d d (4.90)
1
Ni Nj Ni Nj Ni Nj Ni Nj
D12 D33 D11 D33
y x x y y y x x

Estudiando la naturaleza de los distintos términos del integrando se observa que si el

determinante del Jacobiano J es constante, el integrando es un polinomio. Sin
embargo si J no es constante, el integrando es un cociente de polinomios. En el
primer caso la integral puede efectuarse de forma exacta empleando métodos
numéricos adecuados, mientras que en el segundo la integración numérica siempre
es aproximada.

4.9.3 Fuerzas de volumen

Su expresión general es:

Pv NT qv dv (4.91)

Con objeto de simplificar la implementación práctica del método, es habitual

restringir la posible variación de las fuerzas de volumen y limitarla sólo a aquellas que
pueden ser representadas por las propias funciones de interpolación. De esta forma
la variación de las fuerzas de volumen se puede representar mediante una
interpolación de sus valores nodales, empleando las propias funciones de
interpolación:
qv N qev (4.92)

siendo qev los valores nodales de las fuerzas de volumen, que son constantes. Con
esto se obtiene la siguiente expresión del vector de fuerzas nodales equivalentes:

Pv NT N dv qev M qev (4.93)

Considerando la estructura de la matriz N, la matriz M es:

NT1
M NT N dv ... N1 ... Nn t dxdy (4.94)
NTn

con lo que la matriz M se puede dividir en n x n submatrices, que relacionan a los n

nudos entre sí:

M11 ... M1n

M ... ... ... (4.95)
Mn 1 ... Mnn

siendo cada una de ellas:

1 Ni N j 0
ij
M tJ d d (4.96)
1
0 Ni N j

El cálculo de la matriz M no presenta ningún problema, efectuándose en el sistema

local de coordenadas, en el que se conoce la matriz de interpolación N. Todos los
términos que la forman son polinomios, por lo que se puede evaluar de forma exacta
por métodos numéricos.

4.9.4 Fuerzas de superficie

En elasticidad plana, las fuerzas de superficie actúan distribuidas sobre los lados de
los elementos, y pueden definirse de dos formas distintas: como fuerzas distribuidas
por unidad de superficie lateral del elemento qs, o como fuerzas distribuidas por
unidad de longitud qL. La relación entre unas y otras es sencillamente el espesor del
elemento: qL= t qs.
La expresión general del vector de fuerzas nodales equivalentes a las fuerzas de
superficie es:

Ps NT qsds (4.97)

Al igual que con las fuerzas de volumen, se restringe la posible variación de las
fuerzas de superficie sólo a aquellas que pueden ser representadas por las propias
funciones de interpolación. De esta forma las fuerzas de superficie se representan
mediante una interpolación de sus valores nodales, empleando las propias funciones
de interpolación:
qs N qes (4.98)

siendo qes los valores nodales de las fuerzas de superficie, que son constantes. Con
esto queda:

Ps NT N ds qes NT N t dl qes Ms qes (4.99)

donde se ha introducido la matriz Ms que tiene una expresión muy similar a la matriz
M, pero integrando al lado del elemento donde se aplican las fuerzas. Está
compuesta por submatrices del tipo:
1 Ni N j 0
ij
M s t dl (4.100)
1
0 Ni N j

 El valor del diferencial de longitud es: dl 2 dx 2 dy 2

Pero los valores de dx y dy son:
x x
dx d d J 11d J 21d
(4.101)
y y
dy d d J 12d J 22d

 Para una cara donde  es constante será d=0 luego:

dl 2 2
J 21d 2 2
J 22d 2
dl 2
J 21 2
J 22 d (4.102)

 Análogamente, en una cara donde  es constante se obtiene:

dl J 112 J 122 d (4.103)

4.9.5 Elementos triangulares

Todo el desarrollo anterior se ha efectuado considerando como coordenadas locales
las , de los elementos rectangulares, pero todo él es válido para cualquier tipo de
coordenadas locales, y en particular para las coordenadas de área usadas en
elementos de forma triangular.
La mayor parte de las expresiones anteriores son válidas si se cambian
adecuadamente los nombres de las variables locales, aunque hay algunas diferencias.
La primera se refiere al hecho de que al haber tres coordenadas locales de área en el
plano, y sólo dos coordenadas generales  la matriz J sería rectangular; y la
segunda es que los límites de integración ya no son constantes entre -1 y +1, sino que
dependen de las coordenadas.
La forma más sencilla de resolver estos problemas es efectuar un cambio en los
nombres de las variables, en la forma:
L1 L2 (4.104)

La tercera coordenada de área es una variable dependiente:

L3 1 L1 L2 1 (4.105)
Las derivadas de las funciones de interpolación son, por lo tanto:
Ni Ni L1 Ni L2 Ni L3 Ni Ni
(4.106)
L1 L2 L3 L1 L3
Pueden calcularse fácilmente, al conocerse la expresión de Ni en función de Li.
Análogamente se obtiene:
Ni Ni L1 Ni L2 Ni L3 Ni Ni
(4.107)
L1 L2 L3 L2 L3
Utilizando estas ecuaciones, se pueden emplear todas las expresiones ya obtenidas
para los elementos rectangulares.

4.10. CONFORMABILIDAD GEOMÉTRICA

Los elementos isoparamétricos pueden tener lados curvos y un número de nudos
variable en cada uno de ellos, lo que les permite adaptarse bien a la forma del
continuo, pero esta misma cualidad puede dar lugar a problemas de conformabilidad
geométrica.
Huecos entre elementos. Esto puede ocurrir si se unen dos elementos que utilizan
funciones de interpolación diferentes en el contorno de unión, como por ejemplo un
elemento de tres nudos por lado (parabólico) unido a dos elementos lineales (dos
nudos por lado) en un contorno curvo (figura 4.28). La ausencia de huecos en la
discretización se garantiza si los dos elementos contiguos utilizan las mismas
funciones de interpolación en su borde de unión, es decir si emplean el mismo
número de nudos en el lado común.
 Figura 4.28 Huecos e interferencias entre elementos

Deformación excesiva de los elementos. La forma de los elementos viene definida por
la transformación entre las coordenadas locales y las generales, con lo que la forma
de un elemento particular queda determinada por las coordenadas de sus nudos xe.
Unos valores incorrectos para éstas pueden dar lugar a elementos muy deformados,
en los que existirán problemas para efectuar las integrales de área necesarias para
calcular sus propiedades de rigidez.
Ocurren dos efectos simultáneos: algunos puntos interiores del elemento quedan
definidos por dos valores distintos de las coordenadas locales, y por otra parte
algunos puntos interiores en coordenadas locales se salen fuera del contorno del
elemento. En la figura 4.29 se observa un elemento excesivamente deformado en el
que las líneas =0.6 y =0.8 se sale fuera de él, y además se cortan entre sí.

=0.8 =0.6

Figura 4.29 Distorsión excesiva en un elemento curvo

Otro error muy habitual es de la numeración incorrecta de un elemento (figura 4.30):

aunque las coordenadas definan un elemento no deformado, un error en el orden de
los nudos implica una deformación excesiva.
 5 2 1

4
8

3 6 7

Figura 4.30 Numeración incorrecta

Los efectos de distorsión excesiva de los elementos se manifiestan, entre otras

formas, en que la matriz Jacobiana J no es definida positiva en todos los puntos del
elemento. Nótese que esta matriz J debe invertirse para poder calcular las derivadas
de las funciones de interpolación en coordenadas generales y por lo tanto poder
calcular la matriz de rigidez. Por lo tanto la matriz Jacobiana (y su determinante)
monitorizan el estado de distorsión del elemento. De hecho para un elemento
cuadrado el valor de J es constante en todo él e igual a 4 veces su área. A medida
que el elemento se separa de la forma ideal cuadrada el valor de J cambia en algunos
de sus puntos interiores, en la medida que dichos puntos estén situados en zonas
más o menos deformadas respecto a la forma ideal.
Se suelen manejar distintos parámetros para cuantificar la distorsión de los
elementos. Uno de los más sencillos y habituales es:
4J min
DP (4.108)
A
donde Jmin es el valor mínimo del determinante de la Jacobiana en cualquier punto
del elemento y A su área. Normalmente el valor mínimo de J se considera igual al
menor valor obtenido en cualquier punto de integración. Si el elemento es un
paralelogramo el valor de DP=1. Habitualmente se considera que si DP < 0.2 la
distorsión es excesiva.
Puede demostrarse que si las funciones de interpolación son lineales, la condición
necesaria para obtener una transformación correcta es que ningún ángulo interior
sea mayor que 180º. Para funciones de interpolación cuadráticas, además debe
cumplirse que los nudos intermedios se encuentren en el tercio central de la
distancia entre los vértices del lado. Para elementos cúbicos no se puede dar una
regla simple, y hay que comprobar numéricamente el valor del Jacobiano.
4.11. REQUERIMIENTOS DE CONVERGENCIA
La convergencia está garantizada si se cumplen las tres condiciones inicialmente
indicadas de forma general para el MEF, que son, para problemas de elasticidad:
continuidad de tipo C0, y capacidad de representar cualquier estado de tensión
constante.
El requerimiento de continuidad C0 se cumple siempre usando elementos
isoparamétricos con las funciones de interpolación adecuadas, incluso con elementos
muy deformados. Para ello se requiere que se cumpla la misma condición ya indicada
antes para la ausencia de huecos en la discretización, es decir que en dos elementos
adyacentes se usen las mismas funciones de interpolación en el borde de unión entre
ellos, o lo que es lo mismo, que tengan el mismo número de nudos en él.
El requerimiento de representar cualquier estado de tensión constante indica que las
funciones de interpolación deben de permitir representar cualquier valor constante
de la derivada primera del desplazamiento.
Para poder representar dicho valor constante, el campo de desplazamientos debe de
tener como mínimo los términos siguientes:
u 1 2
x y
3
v 4 5
x y
6
(4.109)

ya que adoptando los valores adecuados de i se puede representar cualquier estado

de tensión constante y cualquier movimiento como sólido rígido. Particularizando la
expresión anterior a los nudos se obtiene (en dirección x):
Ui 1
x
2 i
y
3 i
(4.110)
Pero por definición de los elementos ocurre que u N iU i luego, sustituyendo el
valor de Ui
u N iU i Ni 1
N i 2x i N i 3yi (4.111)

Comparando con la expresión inicial se deduce que:

Ni 1 N ixi x N iyi y (4.112)

Las dos últimas expresiones no son otras que las de interpolación de las coordenadas,
pero se añade una nueva, que dice que la suma de las funciones de interpolación
tiene que ser 1, a fin de poder representar cualquier estado de tensión constante.
El desarrollo anterior es para elementos isoparamétricos; para los de tipo
subparamétrico se puede demostrar, de forma análoga, que las funciones N’ usadas
para interpolar las coordenadas, deben de cumplir esta misma relación, con tal de
que las N´ se puedan expresar como una combinación lineal de las funciones N
usadas para interpolar el campo de desplazamientos.
4.12. ELEMENTO CUADRILÁTERO DE CUATRO NUDOS NO CONFORME
El elemento de cuatro nudos isoparamétrico puede tener forma de cuadrilátero, y
admite una variación lineal de los desplazamientos, lo cual le hace muy atractivo para
los problemas de elasticidad plana.
Sin embargo se ha encontrado que aún en ejemplos muy simples, como el de una
viga en voladizo, este elemento se comporta muy mal, y ello es debido a que no es
capaz de representar adecuadamente una deformada curva. En efecto, si se aplican
unos desplazamientos de flexión U a sus nudos extremos (figura 4.31), el elemento se
deforma como se muestra en la parte izquierda de la figura, manteniendo sus lados
rectos, siguiendo la ley:
u U v 0 (4.113)

Sin embargo la forma correcta de deformación de flexión en el plano es la mostrada

en la parte derecha de la figura 4.31, cuya expresión analítica es:
2 2
u U v (1 )LU /(2H ) (1 ) UH /(2L) (4.114)

Esta última ecuación proporciona el valor correcto de la deformación de cortadura

xy=0. Sin embargo la deformación del elemento, dada por la ecuación (4.113),
proporciona un valor de xy =U distinto de 0. Por lo tanto el elemento proporciona
un valor incorrecto de v debido a que en su interior aparece una distorsión de
cortadura xy distinta de cero, llamada cortante parásito. Esto hace que el elemento
sea excesivamente rígido ante una flexión en su plano.

U U U U

xy¹0 xy=0
H

U U

L U U

Figura 4.31 Deformación aproximada y exacta del elemento de 4 nudos

Con objeto de corregir este error Wilson y otros (1973) introdujeron dos nuevos
términos entre las funciones de interpolación del elemento, que representan
precisamente las formas de deformación necesarias para modelizar la deformación
vertical v en el problema de flexión en su plano:
2 2
u N iU i (1 )a1 (1 )a 3 (4.115)
2 2
v NV
i i (1 )a2 (1 )a4
donde ai son unos coeficientes, en principio desconocidos, que se añaden a las
variables nodales del elemento y dan lugar a un vector de grados de libertad
ampliado:
T
e
a
U 1 ... V4 a1 a2 a 3 a 4 (4.116)

Por lo tanto el elemento posee ahora doce grados de libertad en lugar de los ocho
iniciales. La matriz de funciones de interpolación es:
2 2
N 1 ... 0 (1 ) 0 (1 ) 0
N 2 2
(4.117)
0 ... N 4 0 (1 ) 0 (1 )
2
Los coeficientes a1 y a2 representan la colaboración de la forma (1 ) a la
2
deformación total, y los coeficientes a3 y a4 la colaboración de la forma (1 ) . Se
les suele llamar variables no-nodales, al no estar asignadas a ningún nudo en
concreto, como lo están las otras funciones de interpolación (figura 4.32).

v=(1-2) a2 u=(1-2) a3

Figura 4.32 Formas de deformación no nodales

Análogamente, la matriz B tiene el aspecto siguiente:

e
B Ba a (4.118)

Donde Ba es la parte asociada a las variables no nodales a:

Ni
0
x
Ni
Bai 0 i 5 6 (4.119)
y
Ni Ni
y x
Siguiendo el proceso general, la matriz de rigidez Kw que se obtiene es de 12 x 12
pues corresponde a los doce parámetros empleados para representar la deformación
del elemento, y su estructura es la siguiente:
K Kua
KW (4.120)
Kau Kaa
 K BT D B dv Rigidez convencional (8x8) asociada a los e.

Kaa BTa D Ba dv Rigidez añadida (4x4), debida a los parámetros a.

Kua BT D Ba dv Rigidez de acople (8x4) entre los e y los parámetros a.

La nueva matriz de rigidez se puede reducir a una de tamaño 8x8 que considera
únicamente los desplazamientos de los cuatro nudos. Para ello se eliminan los
valores no-nodales ai, mediante una técnica de condensación estática de grados de
libertad.
El elemento así obtenido muestra un comportamiento excelente cuando su forma es
rectangular o de paralelogramo. Sin embargo se comporta mal cuando su forma no
es un paralelogramo y sus lados dejan de ser paralelos. Ello es debido a que las
nuevas funciones de interpolación añadidas lo convierten en un elemento
incompatible: por ejemplo uno cualquiera de los nuevos términos añadidos al campo
de deformaciones puede activarse en un elemento, pero no en su vecino,
incumpliendo el tercer criterio de convergencia.
Para corregir este comportamiento, Taylor y otros propusieron en 1976 un método
para hacer que aun siendo incompatible, el elemento pase el "patch test". La
condición que se impone para ello es que si el elemento está sometido a unos
desplazamientos de los nudos tales que definan un estado de tensión constante, se
cumpla que:
- No se exciten las dos nuevas funciones de interpolación añadidas, es decir que a=0.
- No se produzcan fuerzas sobre los grados de libertad no nodales a.
El elemento así obtenido sigue siendo incompatible, pero se garantiza que es capaz
de representar un estado de tensión constante aun si no es un cuadrilátero, y por
otra parte es capaz asimismo de representar la flexión en su plano. Todo esto hace
que este elemento sea muy utilizado en la práctica, dadas sus excelentes
propiedades y su reducido número de nudos, con lo que genera un número de
ecuaciones relativamente pequeño.

4.13. ELEMENTO PLANO HÍBRIDO DE 4 NUDOS

En un elemento híbrido, se interpolan de manera independiente el campo de
desplazamientos y el campo de tensiones; es decir que este último no se obtiene a
través de la ecuación constitutiva aplicada a las deformaciones unitarias, sino que se
emplea para las tensiones otra interpolación diferente.
El campo de deformaciones en el interior del elemento se aproxima mediante la
expresión habitual, en función de las deformaciones nodales del elemento e :
 e
u N (4.121)
e
Siendo el vector de todas las deformaciones nodales del elemento:
T
e
U 1 V1 U 2 V2 U 3 V3 U 4 V4 (4.122)

y N la matriz de funciones de interpolación:

N1 0 N2 0 N3 0 N4 0
N (4.123)
0 N1 0 N2 0 N3 0 N4

Las funciones de interpolación son bilineales, empleando coordenadas normalizadas:

1
Ni (1 i
)(1 i
) i 1, 4
4
2 1



3 4

Figura 4.33 Elemento híbrido de 4 nudos

Las 4 funciones de interpolación se pueden agrupar en un vector:

N1 1 1 1 1
N2 1 1 1 1 1 1 1 1
(4.124)
N3 4 1 4 1 4 1 4 1
N4 1 1 1 1

Nv a0 a1 a2 a3

Los vectores que contienen las coordenadas x e y de los 4 nudos son:

T T
xN x1 x 2 x 3 x 4 yN y1 y2 y3 y4 (4.125)

La forma del elemento se define mediante la transformación isoparamétrica habitual:

x xTN Nv y yTN Nv (4.126)

La matriz jacobiana de dicha transformación es:

 x y
xTN (a1 a 3 ) yTN (a1 a3 )
J (4.127)
x y xTN (a 2 a 3 ) yTN (a 2 a3 )

Su determinante es lineal, y vale:

J j0 j1 j2
j0 (xTN a1 )(yTN a 2 ) (xTN a 2 )(yTN a1 )
(4.128)
j1 (xTN a1 )(yTN a 3 ) (xTN a 3 )(yTN a1 )
j2 (xTN a 3 )(yTN a 2 ) (xTN a 2 )(yTN a 3 )

4.13.1 Deformaciones unitarias

Las deformaciones unitarias compatibles, producidas por las deformaciones u son:

u
x
x
v
y (4.129)
y
xy
u v
y x
Sustituyendo la interpolación de las deformaciones u se obtiene la expresión habitual
del campo de deformaciones unitarias interpolado, que define la matriz B:
e e
u N B (4.130)

N1 N4
0 ... 0
x x
N1 N4
B N 0 ... 0 (4.131)
y y
N1 N1 N4 N4
...
y x y x
Las derivadas de las funciones de interpolación en coordenadas cartesianas se
obtienen de la forma habitual, mediante la jacobiana de la transformación de
coordenadas:
Ni x y Ni Ni
x x
J (4.132)
Ni x y Ni Ni
y y
4.13.2 Interpolación de tensiones
En la formulación híbrida, el campo de tensiones no se obtiene a partir de las
deformaciones unitarias, sino que se supone para él una cierta variación dentro del
elemento. En este elemento de 4 nudos, el campo de tensiones supuesto
corresponde a un campo lineal en sus coordenadas locales, diferente para cada una
de las tres componentes de la tensión en los ejes locales, en función de 9 parámetros
de ajuste , en la forma:

1 0 0 0 0 0 0 1

0 0 0 1 0 0 0 ... (4.133)
0 0 0 0 0 0 1 9

Este campo de tensiones en coordenadas locales se transforma a las coordenadas

generales por medio de la matriz jacobiana. Esta transformación es:
x xy
(J0 )T J0 (4.134)
yx y

La transformación se efectúa mediante la jacobiana evaluada en el centro del

elemento J0, con objeto de garantizar el cumplimiento del patch test.
x y
J 110 J 120
0
J (4.135)
x y J 210 J 220
0

El campo de tensiones en el sistema cartesiano que se obtiene es:

1 0 0 (J 110 )2 (J 210 )2 (J 110 )2 (J 210 )2 2J 110 J 210 2J 110 J 210
x 1
0 2 0 2 0 2 0 2 0 0 0 0
y 0 1 0 (J ) 12 (J ) 22 (J )12 (J ) 22 2J J12 22 2J J12 22 ... (4.136)
xy 0 0 1 J 110 J 120 J 210 J 220 J 120 J 220 J 210 J 220 AJ AJ 9

Siendo la constante AJ J 110 J 220 J 120 J 210 . Los coeficientes i del ajuste en el sistema
cartesiano son combinación lineal de los coeficientes iniciales i y de los términos de
J0 y son por lo tanto independientes.
El campo de tensiones interpolado debe satisfacer la ecuación de equilibrio, en
sentido débil. Esto se cumple si se impone la condición de que el campo de tensiones
no acumule energía con cualquier estado de deformación unitaria no contenida en la
interpolación de deformaciones unitarias supuesta, denominada deformación
unitaria incompatible in . Dicha condición es:
T
in dA 0 (4.137)
La restricción anterior se puede cumplir eliminando de la expresión de interpolación
de tensiones los términos afectados por los coeficientes 6 a 9, dejando sólo un
término en  y otro en  para cada tensión. De esta forma se obtiene la interpolación
propuesta por Pian y Sumihara que emplea 5 parámetros :

1 0 0 (J 110 )2 (J 210 )2
x 1
0 2 0 2
y 0 1 0 (J ) 12 (J )22 ... (4.138)
xy 0 0 1 J 110 J 120 J 210 J 220 5

En el caso de que no existan deformaciones incompatibles, como es el caso en este

elemento, es posible emplear la interpolación de tensiones inicial completa, con los 9
parámetros, pero los estudios efectuados han demostrado que si se emplean 5
parámetros los resultados obtenidos son más precisos.
Desde el punto de vista de la estabilidad, está demostrado (Pian y Wu, 1988) que, en
un elemento sin deformaciones incompatibles, la condición necesaria para que no
haya modos de energía cero en el elemento es n nu r , siendo n el número de
parámetros empleados para interpolar las tensiones, nu el empleado para interpolar
las deformaciones u y r el número de modos de sólido rígido. Para el caso de un
elemento plano de 4 nudos nu=8 y r=3, con lo que se debe emplear como mínimo
n=5 parámetros.
La ley de interpolación de tensiones en el interior del elemento sufre una última
modificación y se pone en la forma:
1 0 0 J 110 J 110 ( ) J 210 J 210 ( )
x 1
0 0 0 0
y 0 1 0 J J ( 12 12 ) J J (
22 22 ) ... (4.139)
xy 0 0 1 J 110 J 120 ( ) J 210 J 220 ( ) 5

S (4.140)

Siendo S una matriz de tamaño 3 x 5 y  un vector con los 5 parámetros del ajuste.
Las constantes , se introducen en la ley de interpolación con objeto de obtener
matrices desacopladas y representan las coordenadas del centro de gravedad del
elemento:
1 1 1 1
dA Jd d dA Jd d
A A A A
El área del elemento se obtiene empleando la expresión habitual del determinante
de la jacobiana:
1 1
A dA Jd d (j0 j1 j2 )d d 4 j0
1 1
El valor de las constantes es:
1 1 1 1 1 j1
dA Jd d (j0 j1 j2 )d d
A A 1 A 1 3 j0
1 1 1 1 1 j2
dA Jd d (j0 j1 j2 )d d
A A 1 A 1 3 j0

4.13.3 Ecuación de equilibrio

Al ser los campos de deformaciones u y de esfuerzos  independientes, la ecuación
de equilibrio se obtiene empleando un principio variacional, basado en la funcional
de Hellinger-Reissner. Como el campo de desplazamientos u contiene únicamente
desplazamientos compatibles, pues la interpolación de deformaciones cumple las
condiciones de compatibilidad entre elementos, la expresión de dicha funcional es:
1 T T
HR (u, ) D 1
dv uT qv dv uT q S dS (4.141)
v
2 v S

El vector qv define las fuerzas de volumen actuantes sobre el elemento y el qS las

fuerzas de superficie.
La matriz D es la matriz constitutiva del estado plano elástico lineal, cuya inversa en
la expresión anterior proporciona las deformaciones unitarias producidas por el
campo interpolado de tensiones 
La condición estacionaria de la funcional anterior es:
T T T
HR D 1
dv uT qv dv uT qS dS 0
v v S

Sustituyendo las expresiones del campo de deformaciones unitarias  y de las

tensiones interpoladas  se obtiene:
T
HR ST D 1 S eT T
S T
ST e
dv
v
(4.142)
eT T eT T
N q v dv N q S dS 0
v S

Agrupando por una parte los términos correspondientes a los coeficientes  y por
otra los correspondientes a las deformaciones de los nudos queda:

T
HR ST D 1 S dv ST dv e

v v
(4.143)
eT T
Sdv NT q v dv NT q S dS 0
v v S
Definiendo las matrices:

H ST D 1 S dv G ST dv (4.144)
v v

y los vectores de fuerzas nodales equivalentes:

Pv NT qv dv PS NT q S dS (4.145)
v S

la condición estacionaria se puede poner en forma compacta:

T e eT
HR H G GT Pv PS 0 (4.146)

Para que esta condición se cumpla ante cualquier variación de las deformaciones
nodales y de los coeficientes de la interpolación de tensiones , se debe cumplir que:
e
H G 0
(4.147)
GT Pv PS 0

Estas son las ecuaciones de equilibrio del elemento y permiten obtener los valores de
los distintos parámetros e y . De la primera de ellas, se puede obtener
directamente el valor de los coeficientes  de la interpolación de esfuerzos, sin
necesidad de ensamblarlas con los restantes elementos:
e
H 1G (4.148)

Sustituyendo en la segunda ecuación se obtiene la ecuación de equilibrio del

elemento:
GT H 1 G e
Pv PS (4.149)

que define su matriz de rigidez:

Ke GT H 1 G (4.150)

4.14. EJEMPLOS. COMPARACIÓN ENTRE FORMULACIONES

4.14.1 Patch test

Este test ha sido empleado con frecuencia en la bibliografía (Pian y Sumihara, 1984).
En él se analiza un dominio rectangular de dimensiones 10 x 2 y espesor unidad,
sustentado de manera isostática en su cara izquierda. El material tiene E=1500 y
=0.25. Se malla mediante 5 elementos, con un mallado irregular de elementos
distorsionados, como muestra la figura 4.34.
Se consideran dos casos de carga distintos:
 Caso A: un momento puro, producido por dos fuerzas iguales y de sentido
contrario de valor 1000 aplicadas en la cara derecha.
Caso B: carga transversal, formada por dos fuerzas verticales iguales de valor 150
aplicadas en la cara derecha.

2 2 1 1 4 150
1000
B

150
2

A
1000

1 1 2 3 3

Figura 4.34 Modelo para el patch test

El modelo se analiza mediante tres tipos de elementos: el elemento isoparamétrico

en formulación estándar, con 4 y 8 nudos, el elemento incompatible de Wilson y el
elemento híbrido con campo de tensiones impuesto, considerando en éste un
número variable de parámetros para el campo de tensiones. En cada caso se
determinan la deformación vertical del punto A situado en el extremo y la tensión en
dirección X en el punto B junto al apoyo.
La tabla 4.1 muestra los resultados obtenidos y la solución exacta. Se observa que el
elemento estándar de 4 nudos es excesivamente rígido, por su propia naturaleza, y
que la formulación de Wilson mejora la respuesta, aunque todavía está lejos de la
tensión correcta, dada la naturaleza deformada de la malla. En todo caso el elemento
híbrido con 5 parámetros es el de mejores prestaciones para deformación y tensión.
Si se emplean 7 o 9 parámetros, las prestaciones del elemento se degradan y con 9
parámetros, el elemento se comporta de forma similar el bilineal de 4 nudos.

Caso A. Momento Caso B. Carga vertical

Elemento AY xB AY xB

Bilineal estándar de 4 nudos 45.65 -1761 50.95 -2447

Incompatible de Wilson 98.40 -2659 100.31 -3615
Híbrido 5 parámetros 96.18 -3014 98.05 -4073
Híbrido 7 parámetros 92.59 -3041 94.55 -4110
Híbrido 9 parámetros 45.96 -2122 51.27 -2956
Isoparamétrico de 8 nudos 99.98 -2965 103.3 -4842
Exacto 100 -3000 102.6 -4050
Tabla 4.1 Resultados del patch test.
4.14.2 Convergencia
Se analiza un dominio rectangular de dimensiones 150 mm x 50 mm y espesor 10
mm, fijo en su cara izquierda (figura 4.35). El material tiene E=21000 kg/mm2 y =0.3.
Se malla mediante elementos de 4 nudos, con un mallado regular de elementos
cuadrados de tamaño decreciente. El número de elementos en la dimensión menor
se hace variar entre 1 y 7, en la dimensión mayor entre 3 y 21, para un número total
de elementos entre 3 y 147.
Sobre la cara de la derecha se aplica una fuerza vertical de valor total 5000 kg,
aplicada en forma de fuerzas puntuales actuando sobre los distintos nudos existentes
en dicha cara.
Se emplean tres formulaciones distintas: el elemento bilineal estándar de 4 nudos, el
elemento incompatible de Wilson y el elemento híbrido con campo de tensión
impuesto, con ajuste de tensiones de 5 parámetros.

DY 5000
50

X

20
150

Figura 4.35 Ejemplo para verificación de la convergencia

La figura 4.36 muestra el valor de la deformación vertical DY en el borde de la

derecha, donde está aplicada la fuerza, frente al número total de nudos del mallado,
el cual da una idea del tamaño de la malla. Dicha deformación se ha obtenido
promediando las deformaciones de todos los nudos existentes en dicho borde
derecho. Se observa que la formulación estándar produce una solución muy rígida y
de convergencia más lenta que las demás. Los elementos incompatible de Wilson e
híbrido producen la misma solución, por ser una malla no distorsionada, hecho que
ya ha sido demostrado en la bibliografía.
La figura 4.37 muestra el valor de la tensión horizontal en un punto situado en la cara
inferior del dominio, a una distancia de x=20 mm del empotramiento, frente al
número de nudos del mallado. Esta tensión se ha obtenido por promediado simple
de las tensiones en los nudos y posterior interpolación bilineal de los valores
promedio en los nudos. Son aplicables a los resultados las mismas consideraciones
que para la deformación.
 3

Desplazamiento Y
en el extremo
2.8

2.6

2.4
Elemento estándar
Incompatible de Wilson
Elemento híbrido
2.2

Número de nudos
2
50 100 150

Figura 4.36 Convergencia de la deformación al refinar la malla

180

Tensión X promediada
(X=20, Y=0)

160

140

Elemento estándar
Incompatible de Wilson
Elemento híbrido
120

Número de nudos
100
50 100 150

Figura 4.37 Convergencia de la tensión al refinar la malla

4.14.3 Test de Cook
En este ejemplo se analiza un dominio en forma de trapecio irregular cuyas
dimensiones se indican en la figura 4.38. Está fijo en uno de sus lados y sometido a
una carga unidad en el lado opuesto. El material tiene E=1 y =0.33 y el espesor es 1.
Se malla mediante elementos de 4 nudos, con un mallado regular de elementos de
tamaño decreciente, fuertemente distorsionados. Se emplean las tres formulaciones
habituales: el elemento bilineal estándar, el elemento incompatible de Wilson y el
elemento híbrido.

16 P=1

48

48

Figura 4.38. Test de Cook

La tabla 4.2 muestra el valor de la deformación en el centro de la cara donde se

aplica la carga unitaria, para distintos mallados. Puede comprobarse la buena
convergencia de los distintos elementos, incluso para mallados relativamente
groseros.

Mallado Estándar 4 nudos Incompatible de Wilson Híbrido 4 nudos

1x1 6.089 23.688 17.418
2x2 11.775 22.839 21.024
4x4 18.253 23.421 22.968
8x8 22.043 23.763 23.647
16 x 16 23.403 23.883 23.854
24 x 24 23.692 23.915 23.901
32 x 32 23.780 23.929 23.921
Tabla 4.2 Resultados del test de Cook
 5
Elasticidad tridimensional

5.1. INTRODUCCIÓN
Los problemas de elasticidad en tres dimensiones son bastante frecuentes en la
práctica ingenieril, y se presentan sobre todo en elementos que por su proceso de
fabricación, o necesidades funcionales no pueden tener una dimensión mucho menor
que las otras dos. Esto ocurre con piezas fundidas o forjadas (p.e. carcasas de
maquinaria, bancadas de máquinas herramienta, soportes y aparatos de apoyo, etc.),
con elementos estructurales en hormigón para apoyo y cimentación (apoyos de
puentes, cimentaciones de máquinas, obras hidráulicas...), y en general en cualquier
estructura en la que no pueda asumirse la hipótesis de que la tensión en la dirección
del espesor sea nula.
El cálculo de tensiones y deformaciones en un sólido tridimensional es un problema
que no tiene mayor complejidad conceptual que el caso bidimensional, por lo que el
MEF se aplicó desde un principio a este tipo de problemas.
Las ecuaciones diferenciales que controlan el problema son similares a las del
problema bidimensional, pero incluyendo una tercera coordenada z, y un tercer
desplazamiento en dicha dirección w. Estas ecuaciones son de orden m=2. La energía
potencial del sistema contiene a las derivadas primeras de las deformaciones, es decir
n=1, con lo que es suficiente con emplear funciones de interpolación de tipo C 0 para
asegurar la convergencia en este problema.

5.2. CAMPO DE DESPLAZAMIENTOS

Un punto cualquiera del sólido tiene tres desplazamientos u,v,w, que son función de
las coordenadas (x,y,z) del punto, y que se agrupan en un vector:
 u
u v (5.1)
w

Un nudo cualquiera de un elemento tiene tres desplazamientos Ui, Vi, Wi. Todos
ellos se agrupan formando el vector de desplazamientos nodales del elemento:
T
e
U 1 V1 W1 U 2 V2 W2 ... (5.2)

W1
w
v
V1
U1
u

Figura 5.1 Elemento finito tridimensional de 8 nudos

Los desplazamientos se interpolan en función de los desplazamientos nodales

mediante las funciones de interpolación:
u N iU i v NV
i i w NW
i i (5.3)

Esta interpolación se representa en la forma matricial habitual:

e
u N (5.4)

siendo N la matriz de funciones de interpolación, cuya estructura en este caso es:

N1 0 0 ... ... N n 0 0
N 0 N1 0 ... ... 0 Nn 0 (5.5)
0 0 N 1 ... ... 0 0 Nn

5.3. DEFORMACIONES UNITARIAS

En un punto del sólido tridimensional, el vector de deformaciones unitarias contiene
seis términos (Figura 5.2), que corresponden a las tres deformaciones unitarias , y a
las tres deformaciones de cortadura ingenieriles .
Su valor en función de las derivadas de los desplazamientos es, en el caso de
pequeñas deformaciones:
 Z
u
Y X
x
v
x
y
X Y y w
z z
Z (5.6)
xy u v
yz y x
YZ zx
v w
XZ
z y
w u
XY
x z
Figura 5.2 Deformación unitaria

Esta expresión se puede poner en la forma habitual:

0 0
x
0 0
x
y
y
0 0 u
z z
v u (5.7)
xy
0 w
yz y x
zx
0
z y

0
z x
donde se identifica al operador matricial , de tamaño 6x3, que pasa de las
deformaciones u a las deformaciones unitarias. Sustituyendo las deformaciones u en
función de las deformaciones nodales, a través de las funciones de interpolación se
obtiene:
e e
u N B (5.8)

La matriz B relaciona las deformaciones de los nudos con las deformaciones unitarias
en un punto cualquiera del elemento:
B N (5.9)

Dada la estructura de N, la matriz B se puede poner en la forma:

 N1 0 0 ... ... N n 0 0
B N 0 N1 0 ... ... 0 Nn 0 (5.10)
0 0 N 1 ... ... 0 0 Nn

B B1 B2 ... Bn (5.11)

Cada una de las matrices Bi tiene la forma siguiente:

Ni
0 0
x
Ni
0 0
y
Ni 0 0 Ni
0 0
z
Bi 0 Ni 0 (5.12)
Ni Ni
0 0 Ni 0
y x
Ni Ni
0
z y
Ni Ni
0
z x

5.4. RELACIÓN TENSIÓN - DEFORMACIÓN UNITARIA

Las tensiones en un punto cualquiera del dominio están definidas por el tensor de
tensiones en tres dimensiones:
x

z
(5.13)
xy

yz

zx

La ecuación constitutiva del material para un material elástico lineal es:

D( 0
) 0
(5.14)

siendo:
 D la matriz elástica, que para un material elástico lineal es constante y depende
de sólo dos parámetros: el módulo de elasticidad E y el módulo de Poisson . Su
valor es:
 2 0 0 0
2 0 0 0
E
2 0 0 0 (1 )(1 2 )
D (5.15)
0 0 0 0 0 E
0 0 0 0 0 2(1 )
0 0 0 0 0

 0 el vector de las deformaciones unitarias iniciales existentes en el material en el

punto considerado, debidas normalmente a las temperaturas, aunque pueden
incluirse en ellas las debidas a los errores de forma, etc.
 0 son las tensiones iniciales presentes en el material, de valor conocido.

5.5. TIPOS DE ELEMENTOS

El orden de derivación de las funciones de interpolación en la matriz B es la unidad,
por lo que es suficiente con emplear funciones de interpolación con continuidad C0.
Estas funciones se presentan a continuación para los distintos elementos existentes,
empleando sistemas de coordenadas locales adecuados a cada uno de ellos.

5.5.1 Elementos prisma rectangular

Para definir las funciones de interpolación se utilizan tres coordenadas locales , , ,
con límites entre -1 a +1. Cuando se representa en este sistema de ejes local, el
elemento es un cubo de lado 2.




Figura 5.3 Coordenadas locales normalizadas en el espacio

Existe una familia completa de elementos prismáticos de base rectangular, que se

corresponden con los elementos rectangulares planos.
5.5.1.1 Elemento lineal de ocho nudos
Es el correspondiente al rectángulo plano de cuatro nudos. Cada una de sus caras es
un elemento rectangular de cuatro nudos (figura 5.4). Sus funciones de interpolación
son:
1
Ni (1 i
)(1 i
)(1 i
) i 1, 8 (5.16)
8
En esta expresión i,i,i son las coordenadas del nudo i.
3
2
4 



7 6
8

Figura 5.4 Elemento de 8 nudos

5.5.1.2 Elemento cuadrático de veinte nudos

Tiene tres nudos por arista (figura 5.5), por lo que éstas pueden ser parábolas, y cada
una de sus caras tiene ocho nudos, como un elemento plano. Cada cara puede ser
una superficie curva parabólica.
Las funciones de los nudos de esquina son:
1
Ni (1 i
)(1 i
)(1 i
)( i i i
2) (5.17)
8
Las de los nudos intermedios para la cara con i=0 son:
1 2
Ni (1 )(1 i
)(1 i
) (5.18)
4

Figura 5.5 Elemento de 20 nudos

5.5.1.3 Elemento cúbico de 32 nudos
Tiene cuatro nudos por arista y cada cara tiene doce nudos. Sus caras pueden por lo
tanto aproximar superficies de grado 3.

Figura 5.6 Elemento de 32 nudos

5.5.1.4 Elementos lagrangianos

Se pueden desarrollar elementos tridimensionales empleando polinomios de
Lagrange como funciones de interpolación. Las funciones de interpolación se definen
como el producto de tres polinomios de Lagrange, uno en cada dirección.
Los elementos así obtenidos contienen gran cantidad de nudos en las caras y en su
interior, y son muy poco usados en la práctica.

5.5.1.5 Elementos con número de nudos por lado variable

De la misma forma que en el plano, existen en el espacio elementos con número de
nudos variable en sus distintas aristas. Las posibles formas son muy numerosas, y
para generar una función de interpolación se aplica la misma técnica que en el caso
plano, aunque aquí cada nudo de esquina se ve afectado por las tres aristas que
confluyen en él.

5.5.2 Elementos tetraedro

Esta familia de elementos espaciales es la variante tridimensional de los elementos
triangulares planos. Para la definición de las funciones de interpolación se usan unas
coordenadas locales del elemento, que son de volumen (Figura 5.7), y se definen
como:
Volumen(PJKL)
LI (5.19)
Volumen(IJKL)
Existen cuatro coordenadas de este tipo, que evidentemente no son independientes,
ya que se cumple que L1 + L2 + L3 + L4 = 1
En su sistema de ejes local normalizado, los elementos son tetraedros de caras
iguales y volumen unidad.
I

LJ

P
J L

LI
K

Figura 5.7 Coordenadas de volumen

5.5.2.1 Elemento tetraedro de cuatro nudos

Las funciones de interpolación son: N i Li i 1, 4
Este elemento representa en su interior un estado de tensión constante. A pesar de
su pobre representación del campo de tensiones, ha sido muy usado por la facilidad
para generar de forma automática mallas de tetraedros; sin embargo en la actualidad
ha sido desbancado por los elementos hexaedros de cuatro u ocho nudos.
1

2 4

Figura 5.8 Elemento tetraedro lineal de 4 nudos

5.5.2.2 Elemento tetraedro cuadrático de diez nudos

Tiene tres nudos por arista y seis por cara. Las funciones de interpolación de las
esquinas son:
N1 (2L1 1)L1 N2 (2L2 1)L2
(5.20)
N3 (2L3 1)L3 N4 (2L4 1)L4

Para los nudos intermedios son:

N5 4 L1 L2 N6 4 L1 L3 (5.21)
 1

5 7

6
2 4
10

8 9

Figura 5.9 Elemento tetraedro de 10 nudos

5.5.2.3 Elemento tetraedro de veinte nudos

Tiene cuatro nudos por arista y diez por cara. Hay un nudo extra en el centro de cada
cara, al igual que en el elemento plano correspondiente. De esta forma se utiliza un
polinomio completo de orden 3 en las tres coordenadas.

Figura 5.10 Elemento tetraedro de 20 nudos

5.5.3 Elementos prisma triangular

Estos elementos tienen forma de prisma con base triangular. Utilizan un sistema de
coordenadas mixto: para el triángulo de la base se utilizan tres coordenadas de área
L1, L2, L3, y a ellas se añade otra coordenada  para definir la altura respecto de la
base, la cual varía entre -1 y +1 (figura 5.11).
En su sistema local el elemento es un prisma con base triangular de área unidad y
altura igual a 2.
 

Figura 5.11 Coordenadas locales para elementos prisma triangular

Forman esta familia (Figura 5.12) los elementos equivalentes a los ya presentados
antes.
 Prisma triangular lineal de seis nudos, con dos nudos por arista.
 Prisma triangular cuadrático, de 15 nudos en total, con tres nudos por arista.
 Prisma triangular cúbico, de 26 nudos, con cuatro nudos por arista, y un nudo en
el centro de cada una de las dos caras triangulares.

Figura 5.12 Elementos prisma triangular

5.6. FORMULACIÓN ISOPARAMÉTRICA

La formulación isoparamétrica en tres dimensiones sigue exactamente los mismos
pasos que en el caso bidimensional, pero añadiendo la coordenada z.
 La interpolación de coordenadas en el espacio se establece en la forma:
x1
y1
x N1 0 0 N2 0 0 ...
z1
y 0 N1 0 0 N2 0 ... x (5.22)
2
z 0 0 N1 0 0 N2 ... y
2
...

x N xe (5.23)

donde el vector xe agrupa a las coordenadas (x,y,z) de todos los nudos del elemento.
Con esta definición un elemento isoparamétrico espacial tiene una forma geométrica
real que está definida por el tipo de funciones de interpolación que emplee. La forma
real de cada cara o lado está definida por el número de nudos que haya en esa cara o
lado: así los lados con dos nudos son rectos, los lados de tres nudos pueden ser
parábolas y los lados de cuatro nudos pueden ser curvas cúbicas. Las caras con tres
nudos son planas, con cuatro nudos son superficies bilineales, etc.
Es posible utilizar asimismo elementos subparamétricos, que pueden tener cierto
interés desde el punto de vista práctico, ya que pueden emplear unos pocos nudos
(p.e. sólo los cuatro vértices) para definir la forma del elemento, y sin embargo
emplear muchos nudos intermedios en los lados para representar el campo de
deformaciones con gran precisión.
 Las derivadas de las funciones de interpolación respecto de las coordenadas
locales se pueden obtener mediante la regla de derivación en cadena:
Ni x y z Ni Ni
x x
Ni x y z Ni Ni
J (5.24)
y y
Ni x y z Ni Ni
z z
De esta expresión se pueden despejar las derivadas en coordenadas generales:
Ni Ni
x
Ni Ni
J 1
(5.25)
y
Ni Ni
z
El vector de la derecha es conocido sin más que derivar las Ni respecto a 
Conociendo J se pueden obtener de la expresión (5.25) todas las derivadas que
forman la matriz Bi.
 El cálculo de J se efectúa apoyándose en la interpolación de coordenadas:
Ni Ni Ni
xi yi zi

Ni Ni Ni
J xi yi zi (5.26)

Ni Ni Ni
xi yi zi

Esta expresión puede ser evaluada con facilidad, ya que las funciones Ni son
conocidas en función de y xi,yi,zi son las coordenadas de los nudos que definen
la forma del elemento.
 Matriz de rigidez. Su expresión general es:

K BT D B dv (5.27)
v

Todos los términos que aparecen en ella son conocidos. Su evaluación se efectúa de
forma numérica, siendo aplicables las mismas consideraciones que en el caso plano.
El dominio de integración en este caso es un volumen, cuyo elemento diferencial se
expresa como:
dv dx dy dz Jd d d (5.28)

siendo J el determinante de la matriz Jacobiana, que es, en general, una función de

.
 Fuerzas de volumen. Su expresión general es:

Pv NT qv dv (5.29)

Se supone, al igual que en el caso plano, que las fuerzas de volumen se pueden
interpolar a partir de sus valores nodales, empleando las propias funciones de
interpolación:
qv N qev (5.30)

siendo qev los valores nodales de las fuerzas de volumen, que son constantes. Con
esto se obtiene la siguiente expresión del vector de fuerzas nodales equivalentes:

Pv NT N dv qev M qev (5.31)

Considerando la estructura de N, la matriz M es:

 NT1
M NT N dv ... N1 ... Nn dx dy dz (5.32)
NTn

La matriz M se puede dividir en nxn submatrices, que relacionan a los n nudos entre
sí:

M11 ... M1n

M ... ... ... (5.33)
Mn 1 ... Mnn

siendo cada una de ellas:

1
Ni N j 0 0
Mij 0 Ni N j 0 Jd d d (5.34)
1
0 0 Ni N j

El cálculo de la matriz M no presenta ningún problema, efectuándose en el sistema

local de coordenadas, en el que se conoce la matriz de interpolación N. Todos los
términos que la forman son polinomios, por lo que se puede evaluar de forma exacta
por métodos numéricos.
 Fuerzas de superficie. En este caso estas fuerzas actúan sobre una de las caras del
elemento finito, pudiendo tener componentes en las tres direcciones del espacio
(figura 5.13).
La expresión del vector de fuerzas nodales equivalentes a ellas es:

Ps NT qsds (5.35)

Al igual que con las fuerzas de volumen, se supone que las fuerzas de superficie se
representan mediante una interpolación de sus valores nodales, empleando las
propias funciones de interpolación, particularizadas para la cara donde se aplica la
fuerza:
qs N qes (5.36)

siendo qes los valores nodales de las fuerzas de superficie, que son constantes.
 qs =1

dt1 dt2



Figura 5.13 Fuerzas de superficie sobre un elemento espacial

El vector de fuerzas nodales equivalentes es por lo tanto:

Ps NT N ds qse Ms qse (5.37)

donde se ha introducido la matriz Ms que tiene una expresión muy similar a la matriz
M, pero integrando en la cara del elemento donde se aplican las fuerzas. Está
compuesta por submatrices del tipo:

1
Ni N j 0 0
Msij 0 Ni N j 0 ds (5.38)
1
0 0 Ni N j

Por ejemplo, si la fuerza se aplica en la cara =+1 (figura 5.13), el valor del diferencial
de área es: ds dt1 dt2 donde t1 y t2 son dos vectores tangentes a las líneas =Cte
y =Cte, cuyo valor es:

x x
J 11 J 21
y y
dt1 d J 12 d dt2 d J 22 d (5.39)

z J 13 J 23
z

Sustituyendo estos vectores en la expresión de ds se obtiene el valor de éste en

función de las dos coordenadas locales . El dominio de la integración queda así
reducido a un cuadrado de lado 2.
 6
Problemas con simetría de
revolución

6.1. INTRODUCCIÓN
Se dice que un problema tiene simetría de revolución cuando tanto el dominio sólido
del problema, como las cargas aplicadas y las condiciones de ligadura son simétricas
respecto a un eje. La solución del problema, es decir los estados de deformaciones y
tensiones tienen asimismo simetría respecto a dicho eje.
Si se considera una sección del sólido que contiene al eje de revolución, su estado de
deformación queda definido por las dos componentes del desplazamiento u, v
contenidas en dicha sección (figura 6.1). Por lo tanto el problema queda reducido a
un problema en dos dimensiones (radial y axial), aunque el dominio sólido sea
tridimensional.
Por lo tanto el campo de desplazamientos u del problema es:
u(r, z )
u (6.1)
v(r, z )

donde r y z representan respectivamente las coordenadas radial y axial de un punto y

u, v son los correspondientes desplazamientos radial y axial.
Desde el punto de vista de la resolución por el método de los elementos finitos, el
dominio sólido a considerar por cada elemento corresponde a un anillo de revolución
con sección transversal igual a la forma del elemento. Al ser el estado de deformación
bidimensional, se pueden usar las mismas funciones de interpolación utilizadas en el
caso plano para definir el campo de desplazamientos dentro del elemento. Sin
embargo el dominio del elemento sigue siendo el anillo de revolución, por lo que la
energía, la rigidez y las fuerzas nodales equivalentes deberán calcularse para todo el
anillo.
 z z

v
u

Figura 6.1 Problema con simetría de revolución

En un problema de elasticidad plana la energía almacenada está asociada sólo a las

tres componentes de la tensión y la deformación contenidas en el plano del
problema. No contribuyen a la energía ni la tensión ni la deformación
perpendiculares a dicho plano, por ser nulas bien dicha tensión (caso de tensión
plana) o bien la deformación unitaria correspondiente (caso de deformación plana).
Sin embargo, en un problema de revolución todo desplazamiento radial u induce
automáticamente una deformación en dirección circunferencial, y como las tensiones
en esta dirección circunferencial no son nulas, hay que considerar esta cuarta
componente de la tensión y de la deformación unitaria. Este es el punto clave que
marca la diferencia esencial en el tratamiento de los problemas de revolución con
respecto a los problemas planos.

6.2. FUNCIONES DE INTERPOLACIÓN

Al ser la deformación plana, el campo de deformaciones en el interior del elemento
se aproxima mediante la expresión habitual:
u Ni Ui v NV
i i
(6.2)

siendo Ui las deformaciones en dirección radial y Vi las deformaciones en dirección

axial, que forman el vector de deformaciones nodales del elemento:
T
e
U 1 V1 U 2 V2 ... U n Vn (6.3)

donde n es el número de nudos del elemento. La expresión (6.2) puede ponerse en

forma matricial:
e
u N (6.4)
La matriz de funciones de interpolación N tiene dos filas y tantas columnas como
grados de libertad haya entre todos los nudos del elemento. La estructura de esta
matriz es la misma que en un problema de elasticidad plana:
N1 0 N2 0 ... 0 N n 0
N (6.5)
0 N1 0 N2 0 ... 0 Nn

6.3. DEFORMACIONES UNITARIAS

La deformación unitaria en un problema con simetría de revolución tiene cuatro
componentes (figura 6.2), cuya expresión en función de los dos desplazamientos del
punto es:
u
r
r
v
z
z (6.6)
u
rz r
u v
z r
En esta expresión  es la deformación unitaria circunferencial, producida por la
deformación radial, al ser el dominio de revolución. La expresión de  se puede poner
en la forma:
u
0
r r
r
v
0 u
z
z z (6.7)
u 1 v
0
rz r r
u v
z r z r
En esta expresión se identifica al operador matricial  que pasa de las deformaciones
u a las deformaciones unitarias.
 z
r
grz
grz
Z r

t
Y z
f r
X

Figura 6.2 Deformaciones unitarias con simetría de revolución

Sustituyendo las deformaciones u en función de las deformaciones nodales, a través

de las funciones de interpolación se obtiene:
e e
u N B (6.8)

La expresión de la matriz B es:

N1 N2 Nn
0 0 ... 0
r r r
N1 N2 Nn
0 0 ... 0
B N z z z (6.9)
N1 N2 Nn
0 0 ... 0
r r r
N1 N1 N2 N2 Nn Nn
...
z r z r z r
Esta matriz se puede poner en la forma:

B B1 B2 ... Bn (6.10)

siendo cada una de las matrices:

 Ni
0
r
Ni
0
Bi z (6.11)
Ni
0
r
Ni Ni
z r

6.4. ESTADO DE TENSIONES. ECUACIÓN CONSTITUTIVA

El estado de tensiones correspondiente a un problema con simetría de revolución, en
coordenadas cilíndricas, tiene cuatro componentes:
r

z
(6.12)

rz

La ecuación constitutiva, en ausencia de temperaturas, es:

1 0
1 1
r r
1 0
z E (1 ) 1 1 z
(6.13)
(1 )(1 2 ) 1 0
1 1
rz rz
1 2
0 0 0
2(1 )
donde se identifica la matriz elástica D.

6.5. TEMPERATURAS. DEFORMACIONES UNITARIAS INICIALES

La presencia de deformaciones unitarias iniciales introduce un nuevo término en la
ecuación constitutiva:
D( 0
) (6.14)

siendo 0 el vector de las deformaciones unitarias iniciales existentes en el material

en el punto considerado, que son conocidas. Si éstas están producidas por un
incremento de temperatura T, se producen deformaciones unitarias iguales en las
direcciones radial, axial y circunferencial, sin deformación de cortadura.
 T
0Tr

0Tz
T
(6.15)
0
0T T
0Trz 0

6.6. FORMULACIÓN ISOPARAMÉTRICA

La interpolación de coordenadas se emplea de la misma forma que en los problemas
de elasticidad bidimensional, pero aplicada a la definición de las coordenadas radial y
axial:
r1
z1
r N1 0 N2 0 ...
r2 (6.16)
z 0 N1 0 N2 ...
z2
...

x N xe (6.17)

El vector xe agrupa a las coordenadas (r,z) de todos los nudos del elemento.

h x

x N
2

2 r

Figura 6.3 Interpolación de coordenadas

De la misma forma, existen elementos de tipo sub y superparamétrico, en función del

número de nudos usados para definir la forma del elemento. Pueden emplearse
asimismo elementos triangulares, tal y como se indica en el estudio de la elasticidad
bidimensional.
El cálculo de las propiedades del elemento requiere conocer las derivadas de las
funciones de interpolación respecto a las coordenadas generales. Siguiendo las reglas
de la derivación se puede poner la siguiente relación matricial entre las derivadas en
ambos sistemas de coordenadas:
 Ni r z Ni Ni
r J r (6.18)
Ni r z Ni Ni
z z
La matriz J es la matriz Jacobiana de la transformación de coordenadas (r,z) a (xh).
Despejando se obtiene:
Ni Ni
r J 1
(6.19)
Ni Ni
z
El vector de la derecha se puede calcular sin más que derivar las Ni respecto a x,h.
Conociendo J se pueden obtener de la expresión (6.19) todas las derivadas que
forman la matriz Bi.
El cálculo de J se hace apoyándose en la interpolación de coordenadas:
Ni Ni
ri zi
J (6.20)
Ni Ni
ri zi

Esta expresión puede ser evaluada fácilmente, ya que las funciones N son conocidas
en función de xh y ri,zi son las coordenadas de los nudos que definen la forma del
elemento.

6.6.1 Matriz de rigidez

La matriz de rigidez K del elemento es:

K BT D B dv (6.21)

El dominio de integración expresado en las coordenadas locales es:

dv 2 r drdz 2 rJd d (6.22)

El valor del radio se aproxima según la ley de interpolación de coordenadas, con lo

que el diferencial de volumen queda:

dv 2 r drdz 2 N i ri J d d (6.23)

Considerando la estructura de B, la matriz de rigidez K es:

 BT1
K BT D B dv 2 ... D B1 ... Bn J N iri d d (6.24)
BTn

La matriz K se puede dividir en nxn submatrices, que relacionan a los n nudos entre sí:

K11 ... K1n

K ... ... ... (6.25)
Kn 1 ... Knn

Cada una de ellas tiene la expresión:

1
ij
K 2 BTi D B j N iri J d d (6.26)
1

Sobre la naturaleza de estas integrales, pueden hacerse las mismas consideraciones

que en el caso de elasticidad bidimensional. Únicamente debe hacerse notar que el
orden de los polinomios es mayor que en aquel caso, debido a la presencia del
término N iri en el integrando.

6.6.2 Fuerzas de volumen

Las fuerzas de volumen que pueden considerarse en un problema axisimétrico están
distribuidas sobre el volumen de todo el anillo de revolución asociado a cada
elemento finito. Tienen dos componentes, radial y axial, que pueden ser ambas
función de las coordenadas r y z, pero por simetría de revolución no pueden ser
función del ángulo :
qvr
qv q (6.27)
vz

Como ya es habitual, se supone que la variación de las fuerzas de volumen se puede

representar mediante una interpolación de sus valores nodales empleando las
propias funciones de interpolación:
qv N qev (6.28)

siendo qev los valores nodales de las fuerzas de volumen, que son constantes. Con
esto se obtiene la siguiente expresión del vector de fuerzas nodales equivalentes:

Pv NT qv dv NT N dv qev M qev (6.29)

Considerando la estructura de N, la matriz M es:

 NT1
M NT N dv 2 ... N1 ... Nn r dr dz (6.30)
NTn

La matriz M se puede dividir en nxn submatrices, que relacionan a los n nudos entre
sí:

M11 ... M1n

M ... ... ... (6.31)
Mn 1 ... Mnn

Cada una de ellas vale:

1 Ni N j 0
ij
M 2 J N iri d d (6.32)
1
0 Ni N j

El cálculo de la matriz M no presenta ningún problema, efectuándose en el sistema

local de coordenadas, en el que se conoce la matriz de interpolación N. Todos los
términos que la forman son polinomios, por lo que se puede evaluar de forma exacta
de forma numérica.

6.6.3 Fuerzas de superficie

Estas fuerzas corresponden a fuerzas aplicadas sobre la superficie lateral del sólido, y
se transforman en fuerzas aplicadas sobre el lado de los elementos (figura 6.4).
Al igual que con las fuerzas de volumen, las fuerzas de superficie se representan
mediante una interpolación de sus valores nodales, empleando las propias funciones
de interpolación del elemento:
qs N qes (6.33)

siendo qes los valores nodales de las fuerzas de superficie, que son constantes. Con
esto el vector de fuerzas nodales equivalentes resulta ser:

Ps NT N ds qes NT N 2 r dl qes Ms qes (6.34)

donde se ha introducido la matriz Ms que tiene una expresión muy similar a la matriz
M, pero integrando al lado del elemento donde se aplican las fuerzas. Está
compuesta por submatrices del tipo:
Ni N j 0
Msij 2 N iri dl (6.35)
0 Ni N j
El valor del diferencial de longitud se calcula igual que para los elementos
bidimensionales.

qsz

qsr

Figura 6.4 Fuerzas de superficie con simetría de revolución

6.6.4 Fuerzas de línea

Sobre el sólido pueden actuar fuerzas de línea, que deben estar distribuidas sobre
una circunferencia con centro en el eje de revolución (figura 6.5). Para guardar la
simetría de revolución, estas fuerzas deben ser de módulo constante.

z
qLZ 2pRqLZ
R

r
r

Figura 6.5 Fuerzas de línea con simetría de revolución

Sea una fuerza de línea con componentes en las direcciones radial y axial, expresadas
como una fuerza por unidad de longitud circunferencial:
qLr
qL q (6.36)
Lz

Si la fuerza de línea actúa sobre una circunferencia de radio R, la fuerza equivalente

a ella es una fuerza nodal de magnitud:
2 R qLr
PL (6.37)
2 R qLz
Esta fuerza debe aplicarse sobre un nudo de la discretización de elementos finitos
coincidente con la circunferencia de radio R.

6.7. CARGAS SIN SIMETRÍA DE REVOLUCIÓN

Se pretende estudiar ahora un sólido con simetría de revolución, pero sometido a
cargas exteriores que no tienen dicha simetría. En consecuencia las tensiones y las
deformaciones del sólido no tendrán simetría de revolución. Sin embargo quiere
hacerse uso de la simetría del dominio, y emplear elementos finitos planos, en lugar
de elementos tridimensionales.

6.7.1 Campo de desplazamientos

Al no haber simetría de revolución en las cargas, tampoco la hay en las
deformaciones, y ahora el campo de desplazamientos en el interior del sólido es
función de las tres coordenadas r, z, , y tiene tres componentes:
ur
u uz (6.38)
u

Dado que el sólido tiene simetría de revolución, la dependencia de las tres

componentes de u de la coordenada  puede expresarse mediante un desarrollo en
serie de Fourier de la forma:

ur urk1 cos k urk2 sin k

k 0 k 0

uz uzk1 cos k uzk2 sin k (6.39)

k 0 k 0

u u 1k sin k u 2k cos k
k 0 k 0

Los coeficientes del desarrollo se han separado en dos grupos: el grupo con
superíndice 1 corresponde al estado simétrico de deformaciones, y el grupo con
superíndice 2 corresponde al estado antisimétrico. En el estado simétrico los
desplazamientos radial y axial varían según la ley del coseno, que es simétrica
respecto a , mientras que el desplazamiento circunferencial varía según el seno, a fin
de obtener u 0 para =0 y =p, así como el cambio de signo con p<<2p. En el
grupo antisimétrico, las funciones seno y coseno están cambiadas respecto al caso
simétrico.
Todos los coeficientes del desarrollo en serie son sólo funciones de r y , y
constituyen las incógnitas del problema.
6.7.2 Deformaciones unitarias
Las deformaciones unitarias para este problema corresponden a las de la elasticidad
en tres dimensiones, y dada la geometría resulta conveniente manejarlas en
coordenadas cilíndricas:
ur
r
r
uz
z
z
ur 1 u
r r
ur uz
rz (6.40)
z r
u 1 uz
z
z r
1 ur u u
r
r r r
Derivando en las expresiones de los desplazamientos y reagrupando los términos en
seno y coseno, se obtiene la relación entre las deformaciones unitarias y los
coeficientes del desarrollo en serie del campo de desplazamientos. Dicha expresión
se puede poner en forma matricial como:

0 0 0 0 0
r
urk1
r 0 0 0 0 0
z uzk1
z
1 k
0 0 0 0 u 1k
r r
cos k
rz k 0
0 0 0 0 urk2
z r
z
uzk2
k
r 0 0 0 0
r z u 2k
k 1
0 0 0 0
r r r
 0 0 0 0 0
r
urk1
0 0 0 0 0
z uzk1
1 k
0 0 0 0 u 1k
r r
sin k 2
(6.41)
k 0
0 0 0 0 urk
z r
uzk2
k
0 0 0 0
r z u 2k
k 1
0 0 0 0
r r r
C S
k
u12
k
cos k k
u12
k
sin k (6.42)
k 0 k 0

Los operadores Ck y S
k
corresponden a los términos en coseno y seno
respectivamente del desarrollo en serie de . Por su estructura pueden separarse
cada uno de ellos por columnas en dos submatrices, de tamaño 6x3, que
corresponden a los estados 1 y 2, con lo que la ecuación anterior queda:
C1 S1 C2 S2
k
cos k k
sin k u1k k
cos k k
sin k uk2 (6.43)
k 0 k 0

Agrupando entre sí a los operadores de los estados 1 y 2 se obtiene:

0 0 0 0 0
r
r 0 0 0
0 0
z z 0 0 0 urk1
1 k
cos k 0 sin k 0 0 0 uzk1
rz
r r
k 0
k u 1k
z 0 0
z r r z
r
0 0 0 k 1
0
0 0 0 r r r
 0 0 0 0 0
r
0 0 0
0 0
0 0 0 z urk2
1 k
cos k 0 0 0 sin k 0 uzk2 (6.44)
r r
k 0
k u 2k
0 0
r z z r
k 1 0 0 0
0
r r r 0 0 0

1
k
u1k 2
k
uk2 (6.45)
k 0 k 0

1
El primer operador k
corresponde al estado simétrico (estado 1), es de tamaño 6x3,
C S
y está formado por las tres primeras columnas de los operadores k
y k
.
Análogamente el operador k2 , que corresponde al estado antisimétrico (estado 2),
está formado por las tres últimas columnas de los operadores anteriores. Hay que
notar que tanto el operador simétrico como el antisimétrico tienen en su
composición términos en seno y en coseno.

6.7.3 Interpolación de desplazamientos

En las expresiones anteriores se han manejado los desplazamientos u de un punto
cualquiera del sólido. Éstos se interpolan en función de los correspondientes a los
nudos del elemento donde se sitúa el punto, según la expresión habitual:
e
u N (6.46)

donde e contiene a los tres desplazamientos (U en dirección radial, V en dirección

axial y W en dirección circunferencial ) de los n nudos del elemento:
T
e
U 1 V1 W1 U 2 V2 W2 ... (6.47)

y la matriz N es de 3x3n:

N1 0 0 N2 0 0 ...
N 0 N1 0 0 N2 0 ... (6.48)
0 0 N1 0 0 N 2 ...

De la misma forma pueden interpolarse las componentes de cualquiera de los

armónicos que forman la descomposición en serie de la respuesta total, tanto para el
caso simétrico como para el antimétrico:
 e1 e2
u1k N k
uk2 N k
(6.49)

Donde ek1 y ek2 contienen los 3n desplazamientos de los n nudos del elemento e, en
el armónico k, en los casos 1 y 2 respectivamente. Sustituyendo esta expresión en la
de las deformaciones unitarias se obtiene:
1 e1 2 e2
k
N k k
N k
(6.50)
k 0 k 0

e1 e2
B1k k
Bk2 k
(6.51)
k 0 k 0

Esta expresión define las dos matrices B1 y B2 que relacionan los desplazamientos
nodales con las deformaciones unitarias. Hay tantas matrices como armónicos se
emplean en cada uno de los estados 1 y 2. Su expresión es:
c1 s1
B1k k
N cos k k
N sin k Bck1 cos k Bsk1 sin k (6.52)
c2 s2
Bk2 k
N cos k k
N sin k Bck2 cos k Bsk2 sin k (6.53)
Cada una de las matrices B1 y B2 tiene una estructura similar a los operadores
correspondientes, con 3 filas y 3n columnas:
N1
0 0 ... 0 0 0 ...
r
N1 0 0 0 ...
0 0 ...
z 0 0 0 ...
N1 kN 1
Bck1 0 ... Bsk1 0 0 0 ... (6.54)
r r
N1 N1 kN 1 N1
0 ... 0 ...
z r r z
0 0 0 ... kN 1 N1 N1
0 ...
r r r
0 0 0 ...

N1
0 0 0 ... 0 0 ...
r
0 0 0 ... N1
0 0 ...
0 0 0 ... z
N1 kN 1
Bck2 0 0 0 ... Bsk2 0 ... (6.55)
r r
kN 1 N1 N1 N1
0 ... 0 ...
r z z r
kN 1 N1 N1 0 0 0 ...
0 ...
r r r
0 0 0 ...
6.7.4 Estado de tensiones
El vector de tensiones corresponde al estado tridimensional, en coordenadas
cilíndricas:
T

r z rz z r
(6.56)

La ecuación constitutiva tiene la forma D , donde los valores de los coeficientes

de D son los mismos que se emplean en elasticidad tridimensional.

6.7.5 Energía elástica de deformación

Su expresión en régimen lineal es:
1 T 1 T
U dv D dv (6.57)
2 2
Sustituyendo la expresión de las deformaciones unitarias se obtienen tres sumandos
diferentes:
U U1 U2 U 12 (6.58)
siendo:
1 e 1T
U1 p B1pT DBq1dv e1
q (6.59)
2 p q

1 e 2T
U2 p B2pT DBq2dv e2
q (6.60)
2 p q

e 1T
U 12 p B1pT DBq2dv e2
q (6.61)
p q

El primer término corresponde a la energía que se acumula en el caso simétrico 1, el

segundo a la energía en el estado antimétrico 2, y el tercero a la que se acumula
cuando sobre las deformaciones del estado simétrico actúan las tensiones del estado
antimétrico.
A continuación se desarrollan las integrales de los distintos términos.

6.7.5.1 Término 1
Uno cualquiera de los sumandos de estos términos vale:
1 e 1T
1
U pq p Bcp1T DBcq1 cos p cos q dv Bcp1T DBqs 1 cos p sin q dv e1
q
2
1 e 1T
p Bsp1T D Bqc1 sin p cos q dv Bsp1T D Bqs 1 sin p sin q dv e1
q (6.62)
2
Todas las integrales están extendidas al volumen del elemento finito, cuyo elemento
diferencial es dv = r dr dz d.
Dado que las matrices B no dependen de la coordenada circunferencial , es posible
efectuar en primer lugar la integración en dicha coordenada. Esta integración en la
coordenada  puede efectuarse con sencillez si se tiene en cuenta que:
0 para p q 0
sin p sin q d para p q 0 (6.63)
0 para p q

2 para p q 0
cos p cos q d para p q 0 (6.64)
0 para p q

sin p cos q d 0 para todo p, q (6.65)

Con esto la expresión de la energía queda:

1 e 1T
U1 0
2 Bc01T DBc01 dA e1
0
2
(6.66)
1 e 1T c 1T c1 s 1T s1 e1
p B DB dA
p p B DB dA
p p p
2 p 1

siendo dA = r dr dz, y estando ahora todas las integrales extendidas al dominio

plano del elemento.

6.7.5.2 Término 2
Este término tiene una expresión similar a la del término 1, cambiando únicamente el
superíndice 1 por el 2; por lo tanto efectuando el mismo desarrollo se llega a:
1 e 2T T
U2 0
2 Bc02 DBc02 dA e2
0
2
(6.67)
1 e 2T c 2T c2 s 2T s2 e2
p B p DB dAp B p DB dA p p
2 p 1

6.7.5.3 Término 1-2

Uno cualquiera de los sumandos de este término de la energía se puede expandir de
la misma forma que para los términos anteriores, empleando las mismas integrales.
Resulta ser:
1 e 1T
U 12 0
2 Bc01T DBc02 dA e2
0
2
(6.68)
1 e 1T c 1T c2 s 1T s2 e2
p B DB dA
p p B DB dA
p p p
2 p 1
Dada la estructura de las distintas matrices involucradas en los integrandos, resulta
que todas las integrales son nulas, para cualquier valor de p. Por lo tanto este término
de la energía es nulo U12=0.

6.7.6 Matriz de rigidez

Analizando las expresiones de los distintos términos de la energía elástica se
identifican en ellas una serie de matrices de rigidez, que corresponden a los distintos
armónicos en que se ha descompuesto el análisis.
 Caso simétrico 1:
T
K10 2 Bc01 DBc01rdrdz (6.69)
A

T T
K1p Bcp1 DBcp1rdrdz Bsp1 DBsp1rdrdz p 1, (6.70)
A A

 Caso antisimétrico 2:
T
K20 2 Bc02 DBc02rdrdz (6.71)
A

T T
K2p Bcp2 DBcp2rdrdz Bsp2 DBsp2rdrdz p 1, (6.72)
A A

6.7.7 Fuerzas nodales equivalentes

A las fuerzas exteriores actuantes sobre el sólido se les aplica el mismo proceso de
descomposición en serie de Fourier que a los desplazamientos. Para simplificar, en el
desarrollo siguiente se considerarán solamente las fuerzas de volumen, definidas por
sus tres componentes:
qvr
qv qvz (6.73)
qv

Las restantes fuerzas se tratan de un modo similar. La dependencia de las tres

componentes de qv de la coordenada  puede expresarse mediante un desarrollo en
serie de Fourier de la forma:
1 2
qvr qvrk cos k qvrk sin k (6.74)
k 0 k 0

1 2
qvz qvzk cos k qvzk sin k (6.75)
k 0 k 0

qv qv1 k sin k qv2 k cos k (6.76)

k 0 k 0
Estas expresiones pueden ponerse en forma matricial como:

qv A1k q1vk Ak2 qvk

2
(6.77)
k 0 k 0

Las matrices A1 y A2 son:

cos k 0 0 sin k 0 0
A1k 0 cos k 0 Ak2 0 sin k 0 (6.78)
0 0 sin k 0 0 cos k

El campo de desplazamientos puede también expresarse mediante las mismas

matrices:

u A1k u1k Ak2 uk2 (6.79)

k 0 k 0

e introduciendo la ley de interpolación queda:

e1 e2
u A1k N k
Ak2 N k
(6.80)
k 0 k 0

El potencial de las fuerzas de volumen es:

V uT q vdv (6.81)
v

Sustituyendo qv y u se obtienen 4 términos: V V1 V2 V 12 V 21

 Término V1: corresponde al potencial acumulado por la componente simétrica de
las fuerzas sobre la componente simétrica de las deformaciones.
e 1T
V1 p NT A1pT Aq1 q1vqdv (6.82)
p q

Considerando la naturaleza de las matrices A, la integral en  puede efectuarse de

forma independiente, teniendo en cuenta las integrales ya empleadas para el cálculo
de la energía elástica. Se obtiene:
qvr1 0
V1 2 δe01T NT qvz1 0 dA δep1T NT q1vp dA (6.83)
p 1
0

V1 δep1T Pvp1 (6.84)

p 0

donde se han definido los vectores de fuerzas nodales equivalentes

 qvr1 0
Pv10 2 NT qvz1 0 dA Pvp1 NT q1vpdA (6.85)
0

 Término V2: corresponde al potencial acumulado por la componente antisimétrica

de las fuerzas sobre la componente antisimétrica de las deformaciones.
e 2T
V2 p NT A2pT Aq2 qvq
2
dv (6.86)
p q

Considerando la naturaleza de las matrices A, la integral en  puede efectuarse de

forma independiente, teniendo en cuenta las integrales ya empleadas para el cálculo
de la energía elástica. Se obtiene:

0
V2 2 δe02T NT 0 dA δep2T NT q vp
2
dA (6.87)
p 1
qv2 0

V2 δep2T Pvp2 (6.88)

p 0

donde se han definido los vectores de fuerzas nodales equivalentes

0
Pv20 2 NT 0 dA Pvp2 NT q vp
2
dA (6.89)
qv1 0

 Términos V12 y V21. Estos términos corresponden al potencial producido por las
fuerzas del caso simétrico sobre las deformaciones del caso antisimétrico y
viceversa. Dada la naturaleza de las matrices A, estos términos son ambos nulos.

6.7.8 Ecuaciones de equilibrio

El potencial total del elemento es:
e
U V U1 U2 V1 V2 (6.90)
e 1 e 1T e1 1 e 1T e1 1 e 2T e2 1 e 2T e2
0
K10 0 p K1p p 0
K20 0 p K2p p
2 2 p 1 2 2 p 1

δep1T Pvp1 δep2T Pvp2 (6.91)

p 0 p 0

El equilibrio del elemento finito requiere que este potencial sea estacionario para
cualquier variación de sus grados de libertad. De esta condición se obtienen las
ecuaciones de equilibrio del elemento, que son independientes para cada uno de los
armónicos. Las distintas ecuaciones que se obtienen son:
 Caso simétrico 1:
e1
K10 0
P01 (6.92)
e1
K1p p Pp1 p 1, (6.93)

 Caso antisimétrico 2:
e2
K20 0
P02 (6.94)
e2
K2p p Pp2 p 1, (6.95)

En la práctica se emplea un número finito N de armónicos para el análisis, en

cantidad suficiente para representar adecuadamente la respuesta de la estructura. En
consecuencia la ecuación de equilibrio del elemento se expresa mediante 2N
sistemas de ecuaciones, que son independientes unos de otros. Estos sistemas de
ecuaciones se ensamblan con los de los restantes elementos en la forma habitual, y
se resuelven de manera independiente.
Cada uno de los 2N casos a resolver corresponde a un dominio plano situado en el
plano r,z, discretizado con elementos planos, y con tres incógnitas de desplazamiento
en cada nudo Ur, Uz, U. Lógicamente el interés de emplear este método frente a un
mallado tridimensional completo reside en el número de armónicos a emplear para
representar la respuesta del sistema. Si éste es alto, resulta más conveniente estudiar
el problema como tridimensional, pero si es bajo, esta técnica puede resultar más
ventajosa.