CAPÍTULO 31

Método del elemento finito

Hasta aquí hemos empleado métodos por diferencias finitas para resolver ecuaciones
diferenciales parciales. En estos métodos, el dominio de la solución se divide en una
malla con puntos discretos o nodos (figura 31.1b). Entonces, se aplica la EDP en cada
nodo, donde las derivadas parciales se reemplazan por diferencias finitas divididas.
Aunque tal aproximación por puntos es conceptualmente fácil de entender, tiene varias
desventajas. En particular, es difícil de aplicar a sistemas con una geometría irregular,
con condiciones de frontera no usuales o de composición heterogénea.
El método del elemento finito ofrece una alternativa que es más adecuada para tales
sistemas. A diferencia de las técnicas por diferencias finitas, la técnica del elemento
finito divide el dominio de la solución en regiones con formas sencillas o “elementos”
(figura 31.1c). Se puede desarrollar una solución aproximada de la EDP para cada uno
de estos elementos. La solución total se genera uniendo, o “ensamblando”, las soluciones
individuales, teniendo cuidado de asegurar la continuidad de las fronteras entre los
elementos. De este modo, la EDP se satisface por secciones.
Como se observa en la figura 31.1c, el uso de elementos, en lugar de una malla
rectangular, proporciona una mejor aproximación para sistemas con forma irregular.
Además, se pueden generar continuamente valores de las incógnitas a través de todo el
dominio de la solución en lugar de puntos aislados.

FIGURA 31.1
a) Un empaque con geometría irregular y composición no homogénea. b) Un sistema así es muy difícil de modelar con la
técnica por diferencias finitas. Esto se debe al hecho de que se necesitan aproximaciones complicadas en las fronteras del
sistema y en las fronteras entre las regiones de diferente composición. c) Una discretización por elementos finitos es mucho
más adecuada para tales sistemas.

Material A

Material B

Material C

a) b) c)
906 MÉTODO DEL ELEMENTO FINITO

Debido a que una descripción completa va más allá del alcance de este libro, el
presente capítulo ofrece sólo una introducción general al método del elemento finito.
Nuestro objetivo principal es familiarizar al lector con esta técnica y darle a conocer
sus capacidades. Por lo tanto, la siguiente sección ofrece una visión general de los pasos
para la solución de un problema, usando el elemento finito. Después se analizará un
ejemplo sencillo: una barra calentada unidimensional en estado estacionario. Aunque
este ejemplo no usa EDP, nos permite desarrollar y demostrar los principales aspectos
de la técnicas del elemento finito, evitando llegar a factores complicados. Después
podemos analizar algunos problemas con el empleo del método del elemento finito para
resolver EDP.

31.1 EL ENFOQUE GENERAL

Aunque las particularidades varían, la implementación del método del elemento finito
usualmente sigue un procedimiento estándar paso a paso. A continuación se presenta un
panorama general de cada uno de estos pasos. La aplicación de tales pasos a problemas
de ingeniería se desarrollará en las siguientes secciones.

31.1.1 Discretización

Este paso consiste en dividir el dominio de la solución en elementos finitos. En la figu-

ra 31.2 se muestran ejemplos de los elementos empleados en una, dos y tres dimensiones.
Los puntos de intersección de las líneas que forman los lados de los elementos se cono-
cen como nodos, y los mismos lados se denominan líneas o planos nodales.

31.1.2 Ecuaciones de los elementos

El siguiente paso consiste en desarrollar ecuaciones para aproximar la solución de cada

elemento y consta de dos pasos. Primero, se debe elegir una función apropiada con co-
eficientes desconocidos que aproximará la solución. Segundo, se evalúan los coeficien-
tes de modo que la función aproxime la solución de manera óptima.

Elección de las funciones de aproximación. Debido a que son fáciles de manipular

matemáticamente, a menudo se utilizan polinomios para este propósito. En el caso uni-
dimensional, la alternativa más sencilla es un polinomio de primer grado o línea recta.
u(x) = a0 + a1x (31.1)

donde u(x) = la variable dependiente, a 0 y a1 = constantes y x = la variable independien-

te. Esta función debe pasar a través de los valores de u(x) en los puntos extremos del
elemento en x1 y x2. Por lo tanto,
u1 = a0 + a1x1
u2 = a 0 + a1x2
donde u1 = u(x1) y u2 = u(x2). De estas ecuaciones, usando la regla de Cramer, se obtiene
u1 x 2 − u2 x1 u2 − u1
a0 = a1 =
x 2 − x1 x 2 − x1
31.1 EL ENFOQUE GENERAL 907

Elemento línea

a) Unidimensional

Nodo

Línea nodal

Elemento Elemento
cuadrilateral triangular
b) Bidimensional

Elemento
hexaédrico

Plano nodal

c) Tridimensional

FIGURA 31.2
Ejemplos de los elementos empleados en a) una, b) dos y c) tres dimensiones.

Estos resultados se sustituyen en la ecuación (31.1) la cual, después de reagrupar térmi-

nos, se escribe como

u = N1u1 + N2u2 (31.2)

donde

x2 − x
N1 = (31.3)
x 2 − x1
y

x − x1
N2 = (31.4)
x 2 − x1
La ecuación (31.2) se conoce como una función de aproximación, o de forma, y N1 y N2
se denominan funciones de interpolación. Una inspección cuidadosa revela que la ecua-
ción (31.2) es, en realidad, el polinomio de interpolación de primer grado de Lagrange.
Esta ecuación ofrece un medio para predecir valores intermedios (es decir, para inter-
polar) entre valores dados u1 y u2 en los nodos.
908 MÉTODO DEL ELEMENTO FINITO

La figura 31.3 muestra la función de forma junto con las funciones de interpolación
Nodo 1 Nodo 2 correspondientes. Observe que la suma de las funciones de interpolación es igual a uno.
Además, el hecho de que estemos tratando con ecuaciones lineales facilita las ope-
a)
raciones como la diferenciación y la integración. Tales manipulaciones serán importan-
tes en secciones posteriores. La derivada de la ecuación (31.2) es
u1 u
du dN1 dN
u2 = u1 + 2 u2 (31.5)
dx dx dx
De acuerdo con las ecuaciones (31.3) y (31.4), las derivadas de las N se calculan como
b) sigue

1 N1 dN1 1 dN2 1
=− = (31.6)
dx x 2 − x1 dx x 2 − x1

y, por lo tanto, la derivada de u es

c)
du 1
N2 1 = ( −u1 + u2 ) (31.7)
dx x 2 − x1
En otras palabras, es una diferencia dividida que representa la pendiente de la línea
x1 x2 recta que une los nodos.
La integral se expresa como
d)
x2 x2

FIGURA 31.3
∫x1
u dx = ∫x1
N1u1 + N2 u2 dx

b) Una aproximación lineal

Cada uno de los términos del lado derecho es simplemente la integral de un triángulo
o función de forma para
a) un elemento lineal. Las rectángulo con base x2 – x1 y altura u. Es decir,
funciones de interpolación x2 1
correspondientes se
muestran en c) y d).
∫x1
Nu dx =
2
( x 2 − x1 )u

Así, la integral completa es

x2 u1 + u2
∫x1
u dx =
2
( x 2 − x1 ) (31.8)

En otras palabras, esto es simplemente la regla del trapecio.

Obtención de un ajuste óptimo de la función a la solución. Una vez que se ha

elegido la función de interpolación, se debe desarrollar la ecuación que rige el compor-
tamiento del elemento. Esta ecuación representa un ajuste de la función a la solución
de la ecuación diferencial de que se trate. Existen varios métodos para este propósito;
entre los más comunes están el método directo, el método de los residuos ponderados
y el método variacional. Los resultados de todos estos métodos son análogos al ajuste
de curvas. Sin embargo, en lugar de ajustar funciones a datos, estos métodos especifican
relaciones entre las incógnitas de la ecuación (31.2) que satisfacen de manera óptima
la EDP.
31.1 EL ENFOQUE GENERAL 909

Matemáticamente, las ecuaciones del elemento resultante a menudo consisten en un

sistema de ecuaciones algebraicas lineales que puede expresarse en forma matricial,

[k]{u} = {F} (31.9)

donde [k] = una propiedad del elemento o matriz de rigidez, {u} = vector columna de
las incógnitas en los nodos y {F} = vector columna determinado por el efecto de cualquier
influencia externa aplicada a los nodos. Observe que, en algunos casos, las ecuaciones
pueden ser no lineales. Sin embargo, en los ejemplos elementales descritos aquí, así como
en muchos problemas prácticos, los sistemas son lineales.

31.1.3 Ensamble

Una vez obtenidas las ecuaciones de elementos individuales, éstas deben unirse o ensam-
blarse para caracterizar el comportamiento de todo el sistema. El proceso de ensamble
está regido por el concepto de continuidad. Es decir, las soluciones de elementos contiguos
se acoplan, de manera que los valores de las incógnitas (y algunas veces las derivadas)
en sus nodos comunes sean equivalentes. Así, la solución total será continua.
Cuando finalmente todas las versiones individuales de la ecuación (31.9) están
ensambladas, el sistema completo se expresa en forma matricial como

[K]{u′} = {F′} (31.10)

donde [K] = la matriz de propiedades de ensamble y {u′} y {F′} = vectores columna de

las incógnitas y de las fuerzas externas, marcadas con apóstrofos para denotar que son
ensamble de los vectores {u} y {F} de los elementos individuales.

31.1.4 Condiciones de frontera

Antes de resolver la ecuación (31.10) debe modificarse para considerar las condiciones
de frontera del sistema. Dichos ajustes dan como resultado
- -
[k]{u′} = {F ′} (31.11)

donde la barra significa que las condiciones de frontera se han incorporado.

31.1.5 Solución

Las soluciones de la ecuación (31.11) se obtienen con las técnicas que se describieron en
la parte tres, tal como la descomposición LU. En muchos casos, los elementos pueden
configurarse de manera que las ecuaciones resultantes sean bandeadas. Así, es posible
utilizar los esquemas de solución altamente eficientes para estos sistemas.

31.1.6 Procesamiento posterior

Una vez obtenida la solución, ésta se despliega en forma tabular o de manera gráfica.
Además, pueden determinarse las variables secundarias y también mostrarse.
Aunque los pasos anteriores son muy generales, son comunes a la mayoría de las
implementaciones del método del elemento finito. En la siguiente sección ilustraremos
910 MÉTODO DEL ELEMENTO FINITO

cómo se aplican para obtener resultados numéricos de un sistema físico simple (una
barra calentada).

31.2 APLICACIÓN DEL ELEMENTO FINITO EN UNA DIMENSIÓN

En la figura 31.4 se muestra un sistema que puede modelarse mediante la forma unidi-
mensional de la ecuación de Poisson
d 2T
= − f ( x) (31.12)
dx 2
donde f(x) = una función que define una fuente de calor a lo largo de la barra, y donde
los extremos de la barra se mantienen a temperaturas fijas,

T(0, t) = T1

T(L, t) = T2
Observe que ésta no es una ecuación diferencial parcial, sino una EDO con valor en la
frontera. Se usa este modelo sencillo porque nos permitirá introducir el método del elemen-
to finito sin algunas de las complicaciones de una EDP, bidimensional por ejemplo.

EJEMPLO 31.1 Solución analítica para una barra calentada

Planteamiento del problema. Resuelva la ecuación (31.12) para una barra de 10 cm
con las siguientes condiciones de frontera, T(0, t) = 40 y T(10, t) = 200 y una fuente de
calor uniforme de f(x) = 10.

FIGURA 31.4
a) Barra larga y delgada sujeta a condiciones de frontera fijas y una fuente de calor
continua a lo largo de su eje. b) Representación del elemento finito que consta de cuatro
elementos de igual longitud y cinco nodos.

f(x)

T(0, t) T(L, t)
x

x=0 x=L
a)
1 2 3 4 5
1 2 3 4

b)
31.2 APLICACIÓN DEL ELEMENTO FINITO EN UNA DIMENSIÓN 911

Solución. La ecuación a resolver es

d 2T
= −10
dx 2

Suponga una solución de la forma

T = ax2 + bx + c

la cual se deriva dos veces para obtener T″ = 2a. Sustituyendo este resultado en la ecua-
ción diferencial da a = –5. Las condiciones de frontera se utilizan para evaluar los co-
eficientes restantes. Para la primera condición en x = 0,

40 = –5(0)2 + b(0) + c
o c = 40. De manera similar, para la segunda condición,

200 = –5(10)2 + b(10) + 40

de donde se obtiene b = 66. Por lo tanto, la solución final es

T = –5x2 + 66x + 40

Los resultados se grafican en la figura 31.5.

FIGURA 31.5
Distribución de temperatura a lo largo de una placa calentada sujeta a una fuente de calor
uniforme y mantenida a temperaturas fijas en los extremos.

200

100

0
5 10 x
912 MÉTODO DEL ELEMENTO FINITO

31.2.1 Discretización
Nodo 1 Nodo 2
Una configuración simple para modelar el sistema consiste en una serie de elementos de
a) igual longitud (figura 31.4b). Así, el sistema se trata como cuatro elementos de igual
longitud y cinco nodos.
~
T T2
31.2.2 Ecuaciones de los elementos
T1
En la figura 31.6a se muestra un elemento individual. La distribución de temperatura
para el elemento se representa por la función de aproximación
x1 x2
⬃
b) T = N1T1 + N2T2 (31.13)

donde N1 y N2 = funciones de interpolación lineales especificadas por las ecuaciones

(31.3) y (31.4), respectivamente. De esta manera, como se ilustra en la figura 31.6b, la
FIGURA 31.6
a) Un elemento individual.
función de aproximación corresponde a una interpolación lineal entre las dos tempera-
b) Función de aproximación turas nodales.
usada para caracterizar la Como se presentó en la sección 31.1, existen diferentes métodos para desarrollar la
distribución de temperatura ecuación del elemento. En esta sección empleamos dos de ellos. Primero, se usará un
a lo largo del elemento. método directo para el caso sencillo donde f(x) = 0. Posteriormente, debido a su aplica-
bilidad general en ingeniería, dedicamos la mayor parte de la sección al método de los
residuos ponderados.

El método directo. En el caso donde f(x) = 0, se utiliza un método directo para gene-
rar las ecuaciones de los elementos. La relación entre el flujo de calor y el gradiente de
temperatura puede representarse mediante la ley de Fourier:
dT
q = −k ′
dx
donde q = flujo [cal/(cm2 · s)] y k′ = coeficiente de conductividad térmica [cal/(s · cm ·
°C)]. Si se utiliza una función de aproximación lineal para caracterizar la temperatura
del elemento, el flujo de calor hacia el elemento a través del nodo 1 se representa por
T1 − T2
q1 = k ′
x 2 − x1

donde q1 es el flujo de calor en el nodo 1. De manera similar, para el nodo 2,

T2 − T1
q2 = k ′
x 2 − x1

Estas dos ecuaciones expresan la relación de la distribución de la temperatura interna de

los elementos (determinada por las temperaturas nodales) con el flujo de calor en sus
extremos. En consecuencia representan nuestras ecuaciones de los elementos deseadas.
Se simplifican aún más reconociendo que la ley de Fourier se puede utilizar para expre-
sar los flujos de los extremos en términos de los gradientes de temperatura en las fron-
teras. Es decir,
dT ( x1 ) dT ( x 2 )
q1 = − k ′ q2 = k ′
dx dx
31.2 APLICACIÓN DEL ELEMENTO FINITO EN UNA DIMENSIÓN 913

que se sustituyen en las ecuaciones de los elementos para dar

⎧− dT ( x1 ) ⎫
1 ⎡ 1 −1⎤ ⎧T1 ⎫ ⎪⎪ dx ⎪⎪
x 2 − x1 ⎢−1 1 ⎥ ⎨T ⎬ = ⎨ dT ( x ) ⎬ (31.14)
⎣ ⎦ ⎩ 2⎭ ⎪ 2 ⎪
⎪⎩ dx ⎪⎭

Observe que la ecuación (31.14) se presentó en el formato de la ecuación (31.9). Así,

logramos generar una ecuación matricial que describa el comportamiento de un elemen-
to típico de nuestro sistema.
El método directo resulta muy intuitivo. Se utiliza en áreas como la mecánica, para
resolver problemas importantes. Sin embargo, en otros contextos a menudo es difícil, si
no es que imposible, obtener directamente las ecuaciones del elemento finito. En con-
secuencia, como se describe a continuación, se cuenta con técnicas matemáticas más
generales.

El método de los residuos ponderados. La ecuación diferencial (31.12) se reexpre-

sa como
d 2T
0= + f ( x)
dx 2

La solución aproximada [ecuación (31.13)] se sustituye en esta ecuación. Como la ecua-

ción (31.13) no es la solución exacta, el lado izquierdo de la ecuación resultante no será
cero, sino que será igual a un residuo,
d 2 T˜
R = 2 + f ( x) (31.15)
dx
El método de los residuos ponderados (MRP) consiste en encontrar un mínimo para
el residuo, de acuerdo con la fórmula general

∫
D
RWi dD = 0 i = 1, 2,…, m (31.16)

donde D = dominio de la solución y Wi = funciones de ponderación linealmente inde-

pendientes.
Aquí, se tienen múltiples opciones para las funciones de ponderación (cuadro 31.1).
El procedimiento más común para el método del elemento finito consiste en emplear las
funciones de interpolación Ni como las funciones de ponderación. Cuando éstas se sus-
tituyen en la ecuación (31.16), el resultado se conoce como el método de Galerkin,

∫
D
RNi dD = 0 i = 1, 2,…, m

En nuestra barra unidimensional, la ecuación (31.15) se sustituye en esta formulación

para dar

x2 ⎡ d 2 T˜ ⎤
∫
x1
⎢ 2 + f ( x )⎥ Ni dx i = 1, 2
⎢⎣ dx ⎥⎦
914 MÉTODO DEL ELEMENTO FINITO

que se pueden reexpresar como sigue:

x2 d 2 T˜ x2

∫ x1 dx 2
Ni ( x ) dx = − ∫
x1
f ( x ) Ni ( x ) dx i = 1, 2
(31.17)

Ahora, se aplicarán varias manipulaciones matemáticas para simplificar y evaluar

la ecuación (31.17). Una de las más importantes es la simplificación del lado izquierdo
usando la integración por partes. Del cálculo, recuerde que esta operación se expresa
como
b b

∫ u dv = uv a − ∫ v du
b

a a

Cuadro 31.1 Esquemas de residuos alternativos

Se puede elegir entre varias funciones de ponderación para la En el caso de mínimos cuadrados, los coeficientes se ajustan
ecuación (31.16). Cada una representa un procedimiento alter- hasta minimizar la integral del cuadrado del residuo. De manera
nativo para el MRP. que las funciones de ponderación son
En el método de la colocación, elegimos tantas posiciones como
∂R
coeficientes desconocidos existan. Después, se ajustan los coefi- Wi =
cientes hasta que los residuos desaparezcan en cada una de estas ∂ai
posiciones. En consecuencia, la función de aproximación dará
resultados perfectos en las posiciones elegidas, pero en las posi- al sustituir las Wi en la ecuación (31.16), se obtiene
ciones restantes tendremos un residuo diferente de cero. Así, este
∂R
método es parecido a los de interpolación del capítulo 18. Observe
que la colocación corresponde a usar la función de ponderación
∫D
R
∂ai
dD = 0 i = 1, 2,…, n

W = d (x – xi) para i = 1, 2,..., n o

donde n = el número de coeficientes desconocidos y d (x – xi) =
∂
la función delta de Dirac, que es igual a cero en todas partes
excepto en x = xi, donde es igual a 1. ∂ai ∫D
R2 dD = 0 i = 1, 2,…, n

En el método del subdominio, el intervalo se divide en tantos

segmentos, o “subdominios”, como coeficientes desconocidos La comparación de esta formulación con la del capítulo 17
existan. Después, se ajustan los coeficientes hasta que el valor muestra que ésta es la forma continua de la regresión.
promedio del residuo sea cero en cada subdominio. Así, en cada El método de Galerkin emplea las funciones de interpolación
subdominio, la función de ponderación será igual a 1, y la ecua- Ni como funciones de ponderación. Recuerde que estas funciones
ción (31.16) se convierte en siempre suman 1 en cualquier posición en un elemento. En mu-
xi chos problemas, el método de Galerkin da los mismos resultados
∫
x i –1
R dx = 0 para i = 1, 2,…, n que los que se obtienen con los métodos variacionales. En con-
secuencia, ésta es la versión del MRP que se emplea con más
donde xi–1 y xi son las fronteras del subdominio. frecuencia en el análisis del elemento finito.

Si u y v se eligen adecuadamente, la nueva integral en el lado derecho será más fácil de

evaluar que la integral original del lado izquierdo. Esto se puede hacer para el término
del lado izquierdo de la ecuación (31.17), escogiendo Ni (x) como u, y (d2T/dx2) dx como
dv, se obtiene
x2
x2
d 2 T˜ dT˜ x2
dT˜ dNi
∫ x1
Ni ( x ) 2 dx = Ni ( x )
dx dx x1
−
∫x1 dx dx
dx i = 1, 2 (31.18)
31.2 APLICACIÓN DEL ELEMENTO FINITO EN UNA DIMENSIÓN 915

Así, hemos dado el importante paso de bajar el orden en la formulación: de una segunda
a una primera derivada.
A continuación, se evalúa cada uno de los términos que hemos creado en la ecuación
(31.18). Para i = 1, el primer término del lado derecho de la ecuación (31.18) se evalúa
como sigue
x2
dT˜ dT˜ ( x 2 ) dT˜ ( x1 )
N1 ( x ) = N1 ( x 2 ) − N1 ( x1 )
dx x1
dx dx

Sin embargo, de la figura 31.3 recuerde que N1(x2) = 0 y N1(x1) = 1 y, por lo tanto,
x2
dT˜ dT˜ ( x1 )
N1 ( x ) =− (31.19)
dx x1
dx

De manera similar, para i = 2,

x2
dT˜ dT˜ ( x 2 )
N2 ( x ) = (31.20)
dx x1
dx

Así, el primer término en el lado derecho de la ecuación (31.18) representa las condicio-
nes de frontera naturales en los extremos de los elementos.
Ahora, antes de continuar, reagrupemos sustituyendo en la ecuación original los
términos correspondientes por nuestros resultados. Empleamos las ecuaciones (31.18) a
(31.20) para hacer las sustituciones correspondientes en la ecuación (31.17); para i = 1,

x2 dT˜ dN1 dT˜ ( x1 ) x2

∫x1 dx dx
dx = −
dx
+ ∫ x1
f ( x ) N1 ( x ) dx (31.21)

y para i = 2,

x2 dT˜ dN2 dT˜ ( x 2 ) x2

∫x1 dx dx
dx =
dx
+ ∫ x1
f ( x ) N2 ( x ) dx (31.22)

Observe que la integración por partes nos llevó a dos importantes resultados. Pri-
mero, ha incorporado las condiciones de frontera directamente dentro de las ecuaciones
del elemento. Segundo, ha bajado la evaluación de orden superior, de una segunda a una
primera derivada. Este último resultado tiene como consecuencia significativa que las
funciones de aproximación necesitan preservar continuidad de valor, pero no pendiente
en los nodos.
Observe también que ahora podemos comenzar a darles significado físico a cada
uno de los términos que obtuvimos. En el lado derecho de cada ecuación, el primer
término representa una de las condiciones de frontera del elemento; y el segundo es el
efecto de la función de fuerza del sistema, en este caso, la fuente de calor f(x). Como
ahora será evidente, el lado izquierdo representa los mecanismos internos que rigen la
916 MÉTODO DEL ELEMENTO FINITO

distribución de la temperatura del elemento. Es decir, en términos del método del ele-
mento finito, el lado izquierdo será la matriz de propiedad del elemento.
Para ver esto nos concentramos en los términos del lado izquierdo. Para i = 1, el
término es
˜ dN
x2 dT

∫ x1 dx dx
1
dx (31.23)

Recordemos de la sección 31.1.2 que la naturaleza lineal de la función hace que la dife-
renciación y la integración sean sencillas. Si empleamos las ecuaciones (31.6) y (31.7)
para hacer las sustituciones correspondientes en la ecuación (31.23), obtenemos
x1 T1 − T2 1
∫x2 ( x 2 − x1 ) 2
dx =
x1 − x 2
(T1 − T2 ) (31.24)

De manera similar para i = 2 [ecuación (31.22)],

x1 − T1 + T2 1
∫x2 ( x 2 − x1 ) 2
dx =
x1 − x 2
( − T1 + T2 ) (31.25)

Una comparación con la ecuación (31.14) nos muestra que éstas son similares a las
relaciones obtenidas con el método directo usando la ley de Fourier, lo cual se aclara
más al expresar las ecuaciones (31.24) y (31.25) en forma matricial como sigue:

1 ⎡ 1 −1⎤ ⎧T1 ⎫
x 2 − x1 ⎢−1 1 ⎥ ⎨T ⎬
⎣ ⎦ ⎩ 2⎭

Si este resultado se sustituye en las ecuaciones (31.21) y (31.22), y después se expre-

sa en forma matricial, obtenemos la versión final de las ecuaciones de los elementos.

⎧− dT ( x1 ) ⎫ ⎧ 2 f ( x ) N ( x ) dx ⎫
x

1 ⎡ 1 − 1 ⎤ ⎪
⎪ ⎪
dx ⎪ ⎪ x1 ⎪ ∫ 1 ⎪⎪
⎢ ⎥ {T } = ⎨ +
⎬ ⎨ x2 ⎬ (31.26)
x − x1 ⎣−1 1 ⎦ dT ( x 2 ) ⎪ ⎪

2   ⎪
Matriz de rigidez del elemento ⎪⎩ dx ∫
⎪⎭ ⎪⎩ x1



f ( x ) N 2 ( x ) dx
⎭
⎪
⎪
Condición Efectos externos
de frontera

Observe que las ecuaciones del elemento pueden obtenerse no sólo mediante los
métodos directo y de los residuos ponderados, sino también usando el cálculo de varia-
ciones (por ejemplo, véase Allaire, 1985). En el caso presente, este método proporciona
ecuaciones idénticas a las deducidas arriba.

EJEMPLO 31.2 Ecuación del elemento en una barra calentada

Planteamiento del problema. Emplee la ecuación (31.26) para desarrollar las ecua-
ciones del elemento dada una barra de 10 cm, con condiciones en la frontera de T(0, t)
= 40 y T(10, t) = 200 y una fuente de calor uniforme con f(x) = 10. Utilice cuatro ele-
mentos del mismo tamaño con longitud = 2.5 cm.
31.2 APLICACIÓN DEL ELEMENTO FINITO EN UNA DIMENSIÓN 917

Solución. El término de la fuente de calor en el primer renglón de la ecuación (31.26)

se evalúa sustituyendo la ecuación (31.3), e integrando para obtener
2.5 2.5 − x
∫0
10
2.5
dx = 12.5

De manera similar, la ecuación (31.4) se sustituye en el término de la fuente de calor del

segundo renglón de la ecuación (31.26), el cual también se integra para obtener
2.5 x−0
∫0
10
2.5
dx = 12.5

Estos resultados, junto con los valores de los otros parámetros, se emplean para susti-
tuirse en la ecuación (31.26) y así obtener
dT
0.4T1 − 0.4T2 = − ( x1 ) + 12.5
dx

y
dT
–0.4T1 + 0.4T2 = ( x 2 ) + 12.5
dx

31.2.3 Ensamble

Antes de que se ensamblen las ecuaciones del elemento, se debe establecer un esquema
de numeración global que especifique la topología o el arreglo espacial del sistema.
Como en la tabla 31.1, esto define la conectividad de los elementos en la malla. Debido
a que este caso es unidimensional, el esquema de numeración parecerá tan predecible
que resulta trivial. Sin embargo, en problemas para dos y tres dimensiones, tal esquema
es el único medio para especificar qué nodos pertenecen a qué elementos.
Una vez que se ha especificado la topología, la ecuación del elemento (31.26) se
puede escribir para otros elementos usando las coordenadas del sistema. Después, éstas
se agregan (una por una) para ensamblar la matriz de todo el sistema (este proceso se
continuará explorando en la sección 32.4). El proceso se ilustra en la figura 31.7.

TABLA 31.1 Topología del sistema para el esquema de segmentación del elemento
finito de la figura 31.4b.

Número de nodos

Elemento Local Global

1 1 1
2 2
2 1 2
2 3
3 1 3
2 4
4 1 4
2 5
918 MÉTODO DEL ELEMENTO FINITO

⎡ 0.4 –0.4 0 0 0⎤ ⎧T1 ⎫ ⎧–dT (x 1) /dx + 12.5⎫

⎢ ⎥ ⎪ ⎪ ⎪ ⎪
⎢ –0.4 0.4 0 0 0⎥ ⎪T 2 ⎪ ⎪ dT (x 2 ) /dx + 12.5 ⎪
⎪ ⎪ ⎪ ⎪
a) ⎢ 0 0 0 0 0⎥ ⎨0⎬ = ⎨ 0 ⎬
⎢ ⎥ ⎪0⎪ ⎪ ⎪
⎢ 0 0 0 0 0⎥ ⎪ ⎪ ⎪ 0 ⎪
⎢ 0
⎣ 0 0 0 0 ⎥⎦ ⎪⎩ 0 ⎪⎭ ⎪⎩ 0 ⎪⎭

⎡ 0.4 –0.4 +0.4 –0.4 0 0 ⎤ ⎧T1 ⎫ ⎧–dT (x 1) /dx + 12.5⎫

⎢ ⎥⎪ ⎪ ⎪ ⎪
⎢ –0.4 0.4 –0.4 0.4 0 0 ⎥ ⎪T 2 ⎪ ⎪ 12.5 + 12.5 ⎪
⎥ ⎪ ⎪ ⎪ ⎪
b) ⎢ 0 0 0 0 ⎨T 3 ⎬ = ⎨ dT (x 3 ) /dx + 12.5 ⎬
⎢ ⎥⎪ ⎪ ⎪ ⎪
⎢ 0 0 0 0 0 ⎥ ⎪0⎪ ⎪ 0 ⎪
⎢ 0 0 0 0 ⎥
0 ⎦ ⎪⎩ 0 ⎪⎭ ⎪⎩ 0 ⎪⎭
⎣
⎡ 0.4 –0.4 0 0 0 ⎤ ⎧T1 ⎫ ⎧–dT (x 1) /dx + 12.5⎫
⎢ ⎥⎪ ⎪ ⎪ ⎪
⎢ –0.4 0.8 –0.4 0 0 ⎥ ⎪T 2 ⎪ ⎪ 25 ⎪
⎪ ⎪ ⎪ ⎪
c) ⎢ 0 –0.4 0.4 +0.4 –0.4 0 ⎥ ⎨T 3 ⎬ = ⎨ 12.5 + 12.5 ⎬
⎢ ⎥⎪ ⎪ ⎪
⎢ 0 0 0 –0.4 0.4 0 ⎥ ⎪T 4 ⎪ ⎪ dT (x 4 ) /dx + 12.5 ⎪⎪
⎢ 0 0 0 0 0 ⎥⎦ ⎪⎩ 0 ⎪⎭ ⎪⎩ 0 ⎪⎭
⎣
⎡ 0.4 –0.4 0 0 0⎤ ⎧T1 ⎫ ⎧–dT (x 1) /dx + 12.5⎫
⎢ ⎥ ⎪ ⎪ ⎪ ⎪
⎢ –0.4 0.8 0 0 0⎥ ⎪T 2 ⎪ ⎪ 25 ⎪
⎪ ⎪ ⎪ ⎪
d) ⎢ 0 –0.4 –0.4 –0.4 0⎥ ⎨T 3 ⎬ = ⎨ 25 ⎬
⎢ ⎥ ⎪T ⎪ ⎪ 12.5 + 12.5 ⎪
⎢ 0 0 0.8 0.4 +0.4 –0.4 ⎥ ⎪ 4⎪ ⎪ ⎪
⎢ 0
⎣ 0 –0.4 –0.4 0.4 ⎥⎦ ⎪⎩T 5 ⎪⎭ ⎪⎩ dT (x 5 ) /dx + 12.5 ⎪⎭

⎡ 0.4 –0.4 0 0 0⎤ ⎧T1 ⎫ ⎧–dT (x 1) /dx + 12.5⎫

⎢ ⎥ ⎪ ⎪ ⎪ ⎪
⎢ –0.4 0.8 –0.4 0 0⎥ ⎪T 2 ⎪ ⎪ 25 ⎪
FIGURA 31.7 ⎪ ⎪ ⎪ ⎪
e) ⎢ 0 0 0.8 –0.4 0⎥ ⎨T 3 ⎬ = ⎨ 25 ⎬
Ensamble de las ⎢ ⎥ ⎪T ⎪ ⎪ ⎪
ecuaciones de todo ⎢ 0 0 –0.4 0.8 –0.4 ⎥ ⎪ 4⎪
⎪
25
⎪
el sistema.
⎢ 0
⎣ 0 0 –0.4 0.4 ⎥⎦ ⎩⎪T 5 ⎭⎪ ⎪⎩ dT (x 5 ) /dx + 12.5 ⎪⎭

31.2.4 Condiciones en la frontera

Observe que conforme se emsamblan las ecuaciones, se cancelan las condiciones de
frontera internas. Así, el resultado final de {F} en la figura 31.7e tiene condiciones
de frontera sólo para el primero y el último nodos. Ya que T1 y T5 están dados, dichas
condiciones de frontera naturales en los extremos de la barra, dT(x1)/dx y dT(x5)/dx,
representan incógnitas. Por lo tanto, las ecuaciones se reexpresan como sigue:

dT
( x1 ) −0.4T2 = – 3.5
dx
0.8T2 –0.4T3 = 41
–0.4T2 +0.8T3 –0.4T4 = 25 (31.27)
–0.4T3 +0.8T4 = 105
dT
–0.4T4 – ( x5 ) = – 67.5
dx
 31.3 PROBLEMAS BIDIMENSIONALES 919

31.2.5 Solución

De la ecuación (31.27) se obtiene

dT
( x1 ) = 66 T2 = 173.75 T3 = 245
dx
dT
T4 = 253.75 ( x5 ) = −34
dx

31.2.6 Proceso posterior

Los resultados se representan gráficamente. En la figura 31.8 se muestran los resultados

del método del elemento finito, junto con la solución exacta. Observe que el cálculo del
elemento finito capta la tendencia general de la solución exacta y, además, da una coin-
cidencia exacta en los nodos. Sin embargo, existe una discrepancia en el interior de cada
elemento debido a la naturaleza lineal de las funciones de forma.

31.3 PROBLEMAS BIDIMENSIONALES

Aunque la “contabilidad” matemática aumenta de forma notable, la extensión del méto-

do del elemento finito a dos dimensiones es similar, conceptualmente, a los problemas

FIGURA 31.8
Resultados al aplicar el método del elemento finito a una barra calentada. También se
muestra la solución exacta.

T Analítica

Elemento finito
200

100

0
5 10 x
920 MÉTODO DEL ELEMENTO FINITO

unidimensionales analizados hasta ahora. De manera que se siguen los mismos pasos
señalados en la sección 31.1.

31.3.1 Discretización

Comúnmente se emplean elementos sencillos, como triángulos o cuadriláteros, en la

malla del elemento finito para dos dimensiones. En este análisis, nos limitaremos a
elementos triangulares del tipo ilustrado en la figura 31.9.

31.3.2 Ecuaciones del elemento

Tal como en el caso unidimensional, el siguiente paso consiste en desarrollar una ecua-
ción para aproximar la solución del elemento. Para un elemento triangular, la aproxima-
ción más sencilla es el polinomio lineal [compare con la ecuación (31.1)]

u(x, y) = a0 + a1,1x + a1,2y (31.28)

donde u(x, y) = la variable dependiente, las a = coeficientes, x y y = variables indepen-

dientes. Esta función debe pasar a través de los valores de u(x, y) en los nodos del trián-
gulo (x1, y1), (x2, y2) y (x3, y3). Por lo tanto,

u1(x, y) = a0 + a1,1x1 + a1,2y1

u2(x, y) = a0 + a1,1x2 + a1,2y2

u3(x, y) = a0 + a1,1x3 + a1,2y3

o, en forma matricial,

⎡1 x1 y1 ⎤ ⎧ a0 ⎫ ⎧u1 ⎫
⎢1 x ⎪ ⎪ ⎪ ⎪
⎢ 2 y2 ⎥⎥ ⎨ a1,1 ⎬ = ⎨u2 ⎬
⎢⎣1 x3 y3 ⎥⎦ ⎪⎩a1,2 ⎪⎭ ⎪⎩u3 ⎪⎭

FIGURA 31.9
Un elemento triangular.

2
3

1
x
31.3 PROBLEMAS BIDIMENSIONALES 921

de donde se obtiene

1
a0 = [u1 ( x 2 y3 − x3 y2 ) + u2 ( x3 y1 − x1 y3 ) + u3 ( x1 y2 − x 2 y1 )] (31.29)
2 Ae

1
a1,1 = [u1 ( y2 − y3 ) + u2 ( y3 − y1 ) + u3 ( y1 − y2 )] (31.30)
2 Ae

1
a1,2 = [u1 ( x3 − x 2 ) + u2 ( x1 − x3 ) + u3 ( x 2 − x1 )] (31.31)
2 Ae

donde Ae es el área del elemento triangular,

1
Ae = [( x 2 y3 − x3 y2 ) + ( x3 y1 − x1 y3 ) + ( x1 y2 − x 2 y1 )]
2

Las ecuaciones (31.29) a (31.31) se sustituyen en la ecuación (31.28). Después de

reagrupar términos semejantes, el resultado se expresa como sigue:

u = N1u1 + N2u2 + N3u3 (31.32)

donde

1
N1 = [( x 2 y3 − x3 y2 ) + ( y2 − y3 ) x + ( x3 − x 2 ) y]
2 Ae
1
N2 = [( x3 y1 − x1 y3 ) + ( y3 − y1 ) x + ( x1 − x3 ) y]
2 Ae
1
N3 = [( x1 y2 − x 2 y1 ) + ( y1 − y2 ) x + ( x 2 − x1 ) y]
2 Ae

La ecuación (31.32) permite predecir valores intermedios en el elemento, con base

en los valores de sus nodos. En la figura 31.10 se muestra la función de forma junto con
las funciones de interpolación correspondientes. Observe que la suma de las funciones
de interpolación es siempre igual a 1.
Como en el caso unidimensional, hay varios métodos para desarrollar las ecuacio-
nes del elemento, basados en la EDP y en las funciones de aproximación. Las ecuaciones
resultantes son considerablemente más complicadas que la ecuación (31.26). Sin embar-
go, como las funciones de aproximación son normalmente polinomios de grado inferior
como la ecuación (31.28), los términos de la matriz final del elemento consistirán de
polinomios de grado inferior y de constantes.

31.3.3 Condiciones en la frontera y ensamble

La incorporación de condiciones en la frontera y el ensamble de la matriz del sistema

también se hacen un poco más complicados cuando la técnica del elemento finito se
aplique a problemas en dos y tres dimensiones. Sin embargo, como en la deducción de
922 MÉTODO DEL ELEMENTO FINITO

u1

a)
x
u3 u2
y
N1

1
b)
x

0
y 0
N2

c)
0 x

1
0
y
N3

d)
0 x
1

y 0

FIGURA 31.10
a) Una función de aproximación lineal para un elemento triangular. Las funciones
de interpolación correspondientes se muestran en los incisos b) a d).

la matriz del elemento, la dificultad está más relacionada con la mecánica del proceso
que con la complejidad conceptual. Por ejemplo, el establecimiento de la topología del
sistema, que fue trivial para el caso unidimensional, se convierte en un asunto de gran
importancia en los casos de dos y tres dimensiones. En particular, la elección de un
esquema de numeración determinará el bandeado de la matriz del sistema resultante y,
por lo tanto, la eficiencia con la que puede resolverse. En la figura 31.11 se muestra el
esquema desarrollado antes para una placa calentada, y que se resolvió con los métodos
por diferencias finitas en el capítulo 29.
 31.4 RESOLUCIÓN DE EDP CON BIBLIOTECAS Y PAQUETES DE SOFTWARE 923

21 22 23 24 25

26 28 30 32

25 27 29 31
16 17 18 19
20
18 20 22 24

17 19 21 23
11 12 13 14
15
10 12 14 16

9 11 13 15
6 7 8 9
10
2 4 6 8

1 3 5 7

1 2 3 4 5

FIGURA 31.11
El esquema de numeración de los nodos y los elementos de una aproximación por elemento
finito de la placa calentada, que se caracterizó previamente por diferencias finitas en el
capítulo 29.

31.3.4 Solución y procesamiento posterior

Aunque los mecanismos de solución son complicados, la matriz del sistema es tan sólo
un conjunto de n ecuaciones simultáneas que pueden usarse para encontrar los valores
de la variable dependiente en los n nodos. En la figura 31.12 se muestra una solución que
corresponde a la solución por diferencias finitas de la figura 29.5.

31.4 RESOLUCIÓN DE EDP CON BIBLIOTECAS Y PAQUETES

DE SOFTWARE

Bibliotecas y paquetes de software pueden ayudarnos a resolver directamente las EDP.

Sin embargo, como se describe en las siguientes secciones, muchas de las soluciones
están limitadas a problemas sencillos, lo cual es particularmente cierto para los casos
de dos y de tres dimensiones. En tales situaciones, los paquetes genéricos (es decir, los
no desarrollados expresamente para resolver EDP, como los paquetes para el elemento
finito) a menudo están limitados a simples dominios rectangulares.
Aunque esto podría parecer una limitante, los problemas sencillos llegan a tener
gran utilidad desde el punto de vista pedagógico. Ocurre así cuando las herramientas de
visualización de los paquetes se utilizan para desplegar los resultados de los cálculos.

31.4.1 Excel
Aunque Excel no tiene la posibilidad de resolver directamente EDP, es un buen ambien-
te para desarrollar soluciones sencillas para las EDP elípticas. Por ejemplo, la presenta-
924 MÉTODO DEL ELEMENTO FINITO

100 100 100 100 100

75

50
75

50
75

50

0 0 0 0 0

FIGURA 31.12
Distribución de temperatura en una placa calentada, calculada con el método del elemento
finito.

ción ortogonal de las celdas de la hoja de cálculo (figura 31.13b) es análoga a la malla
utilizada en el capítulo 29 para modelar la placa calentada (figura 31.13a).
Como en la figura 31.13b, las condiciones de frontera de Dirichlet pueden introdu-
cirse primero en el contorno del bloque de la celda. La fórmula del método de Liebmann
se implementa al introducir la ecuación (29.11) en una de las celdas del interior (como
la celda B2 de la figura 31.13b). Así, el valor de la celda se calcula en función de las
celdas adyacentes. Luego se copia la celda a las otras celdas interiores. Debido a la na-

FIGURA 31.13
Analogía entre a) una malla 87.5 100 100
rectangular y b) las celdas
de una hoja de cálculo. A B C D E
1 87.5 100 100 100 75
2 75 78.57 76.12 69.64 50
75
3 75 63.17 56.25 52.46 50
4 75 42.86 33.26 33.93 50
5 37.5 0 0 0 25
75
100 + T23 + T12 + 75 B1 + C2 + B3 + A2
T13 = B2 =
4 4
a) Malla b) Hoja de cálculo
 31.4 RESOLUCIÓN DE EDP CON BIBLIOTECAS Y PAQUETES DE SOFTWARE 925

turaleza relativa de la instrucción copiar de Excel, todas las demás celdas serán depen-
dientes de sus celdas adyacentes.
Una vez que usted ha copiado la fórmula, probablemente obtendrá un mensaje de
error: Cannot resolve circular references (No se pueden resolver referencias circulares).
Usted puede rectificar esto yendo al menú de herramientas y seleccionando O(pciones).
Luego seleccione Calcular y verifique el cuadro Iteración. Esto permite que la hoja
de cálculo vuelva a calcular (por omisión, son 100 iteraciones) y seguir el método de
Liebmann iterativamente. Después de esto, presione la tecla F9 para volver a calcular
de forma manual la hoja hasta que las respuestas no varíen, lo cual significa que ha
convergido la solución.
Una vez resuelto el problema, se utilizan las herramientas gráficas de Excel para
visualizar los resultados. En la figura 31.14a se muestra un ejemplo. En tal caso, se tiene
que

• Se usó una malla fina

• Se aisló la frontera inferior
• Se agregó una fuente de calor de 150 a la mitad de la placa (celda E5).

Los resultados numéricos de la figura 31.14a pueden ilustrarse con el asistente para
gráficos de Excel. Las figuras 31.14b y c muestran gráficas de superficies tridimensio-

FIGURA 31.14
a) Solución en Excel de la A B C D E F G H I
ecuación de Poisson para 1 87.5 100.0 100.0 100.0 100.0 100.0 100.0 100.0 75.0
una placa con un extremo 2 75.0 89.2 95.8 99.1 99.7 96.6 89.9 77.6 50.0
inferior aislado y una 3 75.0 86.2 94.7 100.9 103.1 96.7 85.5 70.3 50.0
fuente de calor. b) “Mapa 4 75.0 85.7 96.1 106.7 115.3 101.4 85.2 68.2 50.0
topográfico” y c) ilustración 5 75.0 85.5 97.4 114.3 150.0 108.6 85.6 67.3 50.0
6 75.0 84.0 93.4 103.4 111.6 97.4 81.3 65.6 50.0
tridimensional de las
7 75.0 82.2 88.9 94.2 95.6 88.1 76.6 63.6 50.0
temperaturas.
8 75.0 80.9 85.9 88.9 88.4 82.8 73.5 62.2 50.0
9 75.0 80.4 84.9 87.3 86.3 81.1 72.4 61.7 50.0
a)
S9
S8

S7 160
S6 140

S5 120

S4 100
S9
80 S7
S3
60 S5
S2
40 S3
S1 1 2 3 S1
1 2 3 4 5 6 7 8 9 4 5 6
7 8 9

b) c)
926 MÉTODO DEL ELEMENTO FINITO

nales. Por lo general, la orientación y de éstas es la inversa de la hoja de cálculo. Así, el

extremo superior de las temperaturas más altas (100) normalmente se representará en la
parte inferior de la gráfica. Hemos invertido los valores de y en nuestra hoja antes de
graficar, de modo que las gráficas sean consistentes con la hoja de cálculo.
Advierta cómo las gráficas nos ayudan a visualizar lo que sucederá. El calor fluye
hacia abajo desde la fuente hasta las fronteras, formando una imagen parecida a una
montaña. El calor también fluye hacia abajo desde la frontera con temperatura alta has-
ta los dos extremos laterales. Observe cómo el calor fluye hacia el extremo de baja
temperatura (50). Por último, observe cómo el gradiente de temperatura en la dimensión
y tiende a cero para el extremo inferior aislado (∂T/∂y → 0).

31.4.2 MATLAB

Aunque el paquete MATLAB estándar no tiene grandes capacidades para resolver las
EDP, se pueden desarrollar archivos m y funciones con este propósito. Además, su ca-
pacidad para mostrar imágenes es muy útil, en particular para visualizar problemas
bidimensionales.
Para ilustrar esta capacidad, primero desarrollamos la hoja de cálculo en Excel de
la figura 31.14a. Estos resultados pueden guardarse como un archivo de texto, con un
nombre como plate.txt. Después este archivo se traslada al directorio de MATLAB.
Una vez en MATLAB, se carga este archivo escribiendo
>> load plate.txt

Luego, los gradientes se calculan simplemente así

>> [px,py]=gradient(plate);

Observe que éste es el método más simple para calcular gradientes usando los valores
por omisión de dx = dy = 1. Por lo tanto, serán correctas las direcciones y las magnitu-
des relativas.
Por último, se utilizan una serie de comandos para obtener la gráfica. El comando
contour desarrolla una gráfica de contorno de los datos. El comando clabel agrega
etiquetas de contorno a la gráfica. Finalmente, quiver toma los datos del gradiente y los
añade a la gráfica en forma de flechas.
>> cs=contour(plate);clabel(cs);hold on
>> quiver(-px,-py);hold off

Observe que se ha agregado el signo menos, debido al signo menos de la ley de Fourier
[ecuación (29.4)]. Como se ve en la figura 31.15, la gráfica resultante proporciona una
excelente representación de la solución.
Considere que cualquier archivo que esté en el formato adecuado puede introducir-
se en MATLAB para desplegarse de esta manera. Por ejemplo, el cálculo con IMSL
descrito a continuación podría programarse para generar un archivo que se pueda utili-
zar en MATLAB (o en Excel). Compartir archivos entre herramientas es muy común.
Además, los archivos pueden crearse en un lugar con una herramienta, y transmitirse
vía Internet a otro, donde el archivo pueda usarse con otra herramienta. Éste es uno de
los aspectos interesantes de las aplicaciones numéricas modernas.
31.4 RESOLUCIÓN DE EDP CON BIBLIOTECAS Y PAQUETES DE SOFTWARE 927

100
9 +

8 90+ + 60

7
+
70
6

120 + 80
+ 80 110 +
5 + 140
++
100
4 +

1
1 2 3 4 5 6 7 8 9

FIGURA 31.15
Gráficas de contorno generadas en MATLAB y calculadas con Excel, para la placa
calentada (figura 31.14).

31.4.3 IMSL

IMSL tiene algunas rutinas para resolver EDP (tabla 31.2). En este análisis, nos dedica-
mos a la rutina fps2h. Esta rutina resuelve la ecuación de Poisson o la de Helmholtz en
un rectángulo bidimensional usando una solución rápida de Poisson en una malla uni-
forme.
La subrutina fps2h se implementa con la siguiente instrucción CALL:
CALL FPS2H(PRH,BRH,COEF,NX,NY,AX,BX,AY,BY,IBCT,IORD,U,LDU)

TABLA 31.2 Rutinas IMSL para resolver EDP.

Categoría Rutinas Capacidad

Solución de sistemas de EDP

en una dimensión MOLCH Método de líneas con una base cúbica de Hermite
Solución de una EDP en dos
y tres dimensiones FPS2H Solución de Poisson rápida en dos dimensiones
FPS3H Solución de Poisson rápida en tres dimensiones
928 MÉTODO DEL ELEMENTO FINITO

donde
PRH = FUNCIÓN suministrada por el usuario para evaluar el lado derecho de la
ecuación diferencial parcial. La forma es PRH(X, Y), donde
X = valor de la coordenada X. (Entrada)
Y = valor de la coordenada Y. (Entrada)
PRH debe declararse EXTERNA en el programa de llamada.
BRH = FUNCIÓN suministrada por el usuario para evaluar el lado derecho de las
condiciones de frontera.
La forma es BRHS(ISIDE, X, Y), donde
ISIDE = Número de lado. (Entrada) Véase IBCTY abajo para
la definición de los números laterales.
X = valor de la coordenada X. (Entrada)
Y = valor de la coordenada Y. (Entrada)
BRH debe declararse EXTERNA en el programa de llamada.
COEF = Valor del coeficiente de U en la ecuación diferencial. (Entrada)
NX = Número de líneas de la malla en la dirección X. (Entrada) NX debe ser al
menos 4. Véase el comentario 2 para restricciones adicionales en NX.
NY = Número de líneas de la malla en la dirección Y. (Entrada) NY debe ser al
menos 4. Véase el comentario 2 para restricciones adicionales en NY.
AX = El valor de X a lo largo del lado izquierdo del dominio. (Entrada)
BX = El valor de X a lo largo del lado derecho del dominio. (Entrada)
AY = El valor de Y a lo largo de la parte inferior del dominio. (Entrada)
BY = El valor de Y a lo largo de la parte superior del dominio. (Entrada)
IBCT = Arreglo de tamaño 4 que indica el tipo de condición de frontera en cada
lado del dominio o que la solución es periódica. (Entrada) Los lados están
numerados de 1 a 4 como sigue:
Lado Posición
1—Derecho (X = BX)
2—Inferior (Y = AY)
3—Izquierdo (X = AX)
4—Superior (Y = BY)
Hay tres tipos de condiciones de frontera
IBCTY Condición de frontera
1 El valor de U está dado (Dirichlet)
2 El valor de dU/dX está dado (lados 1 y/o 3). (Neumann)
El valor de dU/dY está dado (lados 2 y/o 4).
3 Periódico.
IORD = Orden de precisión de la aproximación por diferencias finitas. (Entrada)
Puede ser 2 o 4. Normalmente se usa IORD = 4.
U = Arreglo de tamaño NX por NY que contiene la solución en los puntos de la
malla. (Salida)
LDU = Dimensión principal de U exactamente como se especificó en el enunciado
de dimensión del programa de llamada. (Entrada)
 31.4 RESOLUCIÓN DE EDP CON BIBLIOTECAS Y PAQUETES DE SOFTWARE 929

EJEMPLO 31.3 Uso del IMSL para encontrar la temperatura de una placa calentada

Planteamiento del problema. Utilice fps2h para determinar las temperaturas de una
placa cuadrada calentada, con las condiciones de frontera fijas del ejemplo 29.1.

Solución. Un ejemplo de un programa principal en Fortran 90 que usa la función fps2h

para resolver este problema se escribe así:
Program Plate
USE msimsl
IMPLICIT NONE
INTEGER ::ncval, nx, nxtabl, ny, nytabl
PARAMETER (ncval=11, nx=33, nxtabl=5, ny=33, nytabl=5)
INTEGER :: i, ibcty(4), iorder, j, nout
REAL ::QD2VL,ax,ay,brhs,bx,by,coefu,prhs,u(nx,ny),utabl,x,xdata(nx),y,ydata(ny)
EXTERNAL brhs, prhs
ax = 0.0
bx = 40.
ay = 0.0
by = 40.
ibcty(1) = 1
ibcty(2) = 1
lbcty(3) = 1
lbcty(4) = 1
coefu = 0.0
iorder = 4
CALL FPS2H(prhs,brhs,coefu,nx,ny,ax,bx,ay,by,ibcty,iorder,u,nx)
DO i=1, nx
xdata(i) = ax + (bx-ax)*FLOAT(i-1)/FLOAT(nx-1)
END DO
DO j=1, ny
ydata(j) = ay + (by-ay)*FLOAT(j-1)/FLOAT(ny-1)
END DO
CALL UMACH (2, nout)
WRITE (nout,’(8X,A,11X,A,11X,A)’) ‘X’, ‘Y’, ‘U’
DO j=1, nytabl
DO i=1, nxtabl
x = ax + (bx-ax)*FLOAT(i-1)/FLOAT(nxtabl-1)
y = ay + (by-ay)*FLOAT(j-1)/FLOAT(nytabl-1)
utabl = QD2VL(x,y,nx,xdata,ny,ydata,u,nx,.FALSE.)
WRITE (nout.’(4F12.4)’) x, y, utabl
END DO
END DO
END PROGRAM
FUNCTION prhs (x, y)
IMPLICIT NONE
REAL :: prhs, x, y
prhs = 0.0
END FUNCTION
REAL FUNCTION brhs (iside, x, y)
IMPLICIT NONE
INTEGER :: iside
930 MÉTODO DEL ELEMENTO FINITO

REAL :: x, y
IF (iside == 1) then
brhs = 50.
ELSEIF (iside == 2) THEN
brhs = 0.
ELSEIF (iside == 3) THEN
brhs = 75.
ELSE
brhs = 100.
END IF
END FUNCTION

Una corrida de ejemplo proporciona la siguiente salida:

X y u x y u
.0000 .0000 37.5000 30.0000 20.0000 52.3849
10.0000 .0000 .0000 40.0000 20.0000 50.0000
20.0000 .0000 .0000 .0000 30.0000 75.0000
30.0000 .0000 .0000 10.0000 30.0000 79.0032
40.0000 .0000 25.0000 20.0000 30.0000 76.8058
.0000 10.0000 75.0000 30.0000 30.0000 69.9017
10.0000 10.0000 42.5976 40.0000 30.0000 50.0000
20.0000 10.0000 32.2945 .0000 40.0000 87.5000
30.0000 10.0000 33.4962 10.0000 40.0000 100.0000
40.0000 10.0000 50.0000 20.0000 40.0000 100.0000
.0000 20.0000 75.0000 30.0000 40.0000 100.0000
10.0000 20.0000 63.5128 40.0000 40.0000 75.0000
20.0000 20.0000 56.2493

PROBLEMAS

31.1 Repita el ejemplo 31.1, pero para T(0, t) = 75 y T(10, t) = donde Ac = área de la sección transversal, E = módulo de Young,
150, y una fuente uniforme de calor de 15. u = deflexión, y x = distancia medida a lo largo de la longitud de
31.2 Repita el ejemplo 31.2, pero para condiciones de frontera la barra. Si la barra está fija rígidamente (u = 0) por ambos ex-
de T(0, t) = 75 y T(10, t) = 150, y una fuente de calor de 15. tremos, use el método del elemento finito para modelar sus de-
31.3 Aplique los resultados del problema 31.2 para calcular la flexiones para Ac = 0.1 m2, E = 200 × 109 N/m2, L = 10 m, y P(x)
distribución de temperatura para la barra completa, con el uso = 1 000 N/m. Emplee un valor de ∆x = 2 m.
del enfoque del elemento finito. 31.6 Desarrolle un programa amigable para el usuario a fin de
31.4 Utilice el método de Galerkin para desarrollar una modelar la distribución de estado estable de la temperatura en
ecuación de elemento para una versión de estado estable de una barra con fuente de calor constante, con el método del ele-
la ecuación de advección-difusión descrita en el problema mento finito. Elabore el programa de modo que se utilicen nodos
30.7. Exprese el resultado final en el formato de la ecuación irregularmente espaciados.
(31.26) de modo que cada término tenga una interpretación fí- 31.7 Utilice Excel para realizar el mismo cálculo que en la fi-
sica. gura 31.14, pero aísle el borde del lado derecho y agregue una
31.5 El modelo siguiente es una versión de la ecuación de Pois- fuente de calor de –150 en la celda C7.
son que ocurre en la mecánica para la deflexión vertical de una 31.8 Emplee MATLAB para desarrollar una gráfica de contorno
barra con una carga distribuida P(x): con flechas de flujo para la solución en Excel del problema 31.7.
31.9 Use Excel para modelar la distribución de temperatura de
∂ 2u la placa que se muestra en la figura P31.9. La placa tiene 0.02
Ac E = P( x ) m de espesor y una conductividad térmica de 3 W/(m · °C).
∂x 2
 PROBLEMAS 931

en 5 elementos (6 nodos de 10 cm de largo cada uno). El extremo

100⬚C izquierdo de la barra tiene un gradiente fijo de temperatura y la
temperatura es variable. El extremo derecho tiene una tempera-
2m
–100 W/m2 tura fija y el gradiente es variable. La fuente de calor f(x) tiene
un valor constante. Así, las condiciones son
75⬚C 1m 0.4 m 25⬚C
dT
= 0.25° C/m T x = 50 = 100° C f ( x ) = 30 W/cm
∂x x =0
0.6 m
Desarrolle las ecuaciones nodales que deben resolverse para las
50⬚C temperaturas y los gradientes de temperatura en cada uno de los
seis nodos. Acople las ecuaciones, inserte las condiciones de
frontera y resuelva el conjunto resultante para las incógnitas.
Figura P31.9 31.13 Encuentre la distribución de temperatura en una barra
(véase la figura P31.13) con generación interna de calor, con el
método del elemento finito. Obtenga las ecuaciones de elemen-
tos nodales con el uso de la conducción de calor de Fourier.
31.10 Use MATLAB para desarrollar una gráfica de contorno
dT
con flechas de flujo para la solución con Excel del problema qk = − kA
31.9. ∂x
31.11 Emplee IMSL para llevar a cabo el mismo cálculo que en y las relaciones de conservación de calor
el ejemplo 31.3, pero aísle el borde inferior de la placa.
31.12 Encuentre la distribución de temperatura en una barra
(véase la figura P31.12) con generación de calor interno, por ∑ [q + f ( x)] = 0
k

medio del método del elemento finito. Obtenga las ecuaciones

nodales de elemento con el uso de la conducción de calor de donde qk = flujo de calor (W), k = conductividad térmica
Fourier. (W/(m · °C)), A = área de la sección transversal (m2) y f(x) =
dT fuente de calor (W/cm). La barra mide 50 cm de longitud, la
qk = − kA coordenada x es cero en el extremo izquierdo y positiva hacia la
∂x
derecha. La barra también tiene bloqueos lineales en x = 0 y en
y las relaciones de conservación del calor
x = 50, con un valor de kA = 100 y 50 W m/°C, respectivamente.
∑ [q + f ( x)] = 0
k
Divida la barra en 5 elementos (6 nodos, cada uno de 10 cm de
largo). Ambos extremos de la barra tienen temperaturas fijas. La
donde qk = flujo de calor (W), k = conductividad térmica fuente de calor f(x) tiene un valor constante. Así, las condiciones
(W/(m · °C)), A = área de la sección transversal (m2), y f(x) = son
fuente de calor (W/cm). La barra tiene un valor de kA = 100 W
m/°C. La barra mide 50 cm de largo, en el extremo izquierdo la T x =0
= 100°C T x = 50
= 50°C f ( x ) = 30 W/cm
coordenada x es cero y positiva hacia la derecha. Divida la barra

Figura P31.13

Figura P31.12
kA = 100 W·m/⬚C
f(x) = 30 W/cm
kA = 50 W·m/⬚C
kA = 100 W/m ·⬚C
dT f(x) = 30 W/cm x
–– = 0.25⬚C/m
dx x=0
Tx=50 = 50⬚C
50 cm
x T x=50 = 100⬚C 50 cm
Tx=0 = 100⬚C
932 MÉTODO DEL ELEMENTO FINITO

Las áreas bloqueadas deben tratarse como si fueran constantes nodales que deben resolverse para las temperaturas y los gra-
en la longitud de un elemento. Por tanto, promedie los valores dientes de temperatura en cada uno de los seis nodos. Ensamble
kA en cada extremo del nodo y tome el promedio como una las ecuaciones, inserte las condiciones de frontera y resuelva
constante en el nodo. Desarrolle las ecuaciones de elementos para el conjunto resultante de incógnitas.