MEF035A165

2.
Modelación y Simulación Numérica de Problemas de Conducción de

Calor.
En la Ingeniería Mecánica, así como en otras ramas de la ciencia e ingeniería se tiene como propósito
esencial la interpretación de fenómenos reales observados, de forma que se puedan crear modelos que
los expliquen y permitan hacer predicciones, y con lo anterior se pretende modificar o controlar su
comportamiento, lo cual es de suma importancia.
Una vez que se puede discriminar de la fenomenología de un problema específico, los elementos
que influyen en el desarrollo de cierto comportamiento del fenómeno, por ejemplo: la fluencia del
material de un elemento de una máquina, la turbulencia en el flujo de fluidos en tuberias, etc.,
se procede a la elaboración de un modelo. El modelo creado refleja la realidad física, dentro de
ciertos límites, en general los modelos por sus limitaciones son “aproximaciones” de la realidad. Los
tipos de modelos de interés son aquellos que nos pueden permitir estudiar el fenómeno, y mediante
su simulación o experimentación permiten conocer la evolución temporal y la variación espacial de
aquellas magnitudes físicas que caracterizan al fenómeno.
Los modelos que inicialmente se considerán se denominan: modelos físicos y modelos experimentales.
Los modelos físicos son aquellos que permiten estudiar los fenómenos mediante “simulación numérica”
(mecánica computacional), y los modelos experimentales permiten estudiar los fenómenos mediante
“experimentación física” (mecánica experimental).
Para un fenómeno específico tenemos entre sus modelos físicos y experimentales tenemos una in-
terrelación dada por ser representaciones del mismo fenómeno, por lo tanto ambos modelos deben
intercambiar información entre sí, de manera que los modelos sean perfectibles, dándose así, una
aproximación significativa al fenómeno en estudio.
La resolución de un modelo físico no es inmediata ni sencilla. Técnicas matemáticas de alto nivel son
necesarias para realizar análisis cualitativos y cuantitativos del problema [Ma 1985], [Ba-Pr 1986],
[Ro 1970], [Na-Se 1982], [Br 1984], etc., así como para generar diversos modelos numéricos [Al 1997],
[Gl 1984], [Gl-Li-Tr 1981], [Fo-Gl 1983], etc., implantables en computadora, por medio de los cuales
se obtendrán finalmente soluciones numéricas discretas del problema en consideración, dentro de
ciertos rangos de precisión e incertidumbre.
La metodología actual para resolver un sistema de ecuaciones de campo (que forman el modelo
físico), centra su atención en el desarrollo de principios variacionales equivalentes al mismo, desde
el punto de vista de las soluciones, cuya estructura intrínseca ha probado ser el marco natural para
los métodos de discretización. Bajo estos últimos se entiende el proceso de replantear en forma
aproximada el problema dentro de espacios de dimensión finita, logrando así sistemas algebraicos
que modelan discretamente el fenómeno.
Existen fenómenos que directamente se caracterizan por poseer un número finito de grados de libertad,
lo cual nos conduce a modelos discretos, en forma natural.
Una vez modelado discretamente el fenómeno, la selección de métodos numéricos apropiados para su
resolución tiene lugar, siendo el siguiente paso su implantación en computadora.
Para el modelo experimental que representa un fenómeno real, se requiere obtener los modelos anál-
ogos del fenómeno, los cuales son interpretaciones del fenómeno en estudio. Estos modelos son
caracterizados por las propiedades de fenómenos físicos conocidos, como lo son: el comportamiento
35
óptico de la luz, la corriente eléctrica, propiedades mecánicas de algunos materiales, etc.
De acuerdo a la analogía o analogias adecuadas para el fenómeno en estudio, tenemos, la creación
de los modelos de ensayo o prueba. De los modelos de ensayo tenemos en la experimentación la
posibilidad de hacer el análisis cualitativo del fenómeno, vgr. la existencia de solución, si existe
periodicidad del fenómeno, etc.
De los modelos de ensayo creados, pueden representar fenómenos continuos o discretos. En el caso
continuo para realizar una cuantificación se determinan puntos de interés de estudio, los cuales nos
servirán para hacer discreto el modelo de ensayo, y así poder cuantificar los parámetros o variables,
asequibles, que son derivados de los modelos análogos mismos.
De los datos obtenidos en la discretización del modelo de ensayo, por lo general requieren un proce-
samiento computacional adecuado para dar la interpretación del fenómeno en estudio.
El tratamiento de interpretación del comportamiento de un fenómeno tanto mediante modelos físicos
como modelos experimentales en forma conjunta, hace que se tenga una visión unificada de una
estructura conceptual y metodológica de los problemas propios de la ingeniería mecánica.
Lo antes expuesto, es resumido en la figura 1.
Figura 1. Modelación de Problemas.

El presente texto trata de la simulación numérica de problemas de la Ingeniería Mecánica.
Nuestro objetivo general es estudiar simulaciones numéricas de problemas de Ingeniería Mecánica
mediante técnicas de elemento finito.
Para ello se estudiarán métodos que permitan realizar la modelación de manera sistemática. Los
problemas de interés son tanto lineales como no lineales, en una o varias dimensiones espaciales y en
general problemas no potenciales.
36
De la figura 1, obtenemos el esquema conceptual de la simulación numérica de problemas de la
Ingeniería Mecánica, asociando las fases de la modelación con las áreas de conocimiento necesarias
para el estudio e implantación sistemática de soluciones mediante la simulación computacional.
Figura 2. Esquema Conceptual.
En el esquema conceptual observamos los conocimientos requeridos, iniciando con la Mecánica del
Medio Continuo [Gu 1981] y la Termodinámica [Ow 1984], mediante estas teorías, podemos expresar
los modelos físicos correspondientes a diversos problemas típicos de la ingeniería [Du-Li 1972], la
generación de los modelos matemáticos correspondientes a los problemas de valores iniciales y de
frontera que representan el modelo físico, pueden ser obtenidos a través de diversas metodologias.
En el caso de problemas potenciales [Ek-Te 1976], [Gl 1984], [Gl-Li-Tr 1981], [Fo-Gl 1983], [Al 1987],
se tiene la teoría de dualidad-perturbación, dualidad-penalización, etc., y para problemas no po-
tenciales se tiene la metodología de subdiferenciales [Al 1989], que junto con la metodología de
resolventes [Al 1992], [Al 1997] generalizan las metodologias para problemas potenciales.
Para la aproximación del problema continuo se hace uso de la aproximación de los espacios de solución
mediante las aproximaciones internas o externas [Te 1973].
La concreción de los modelos discretos, se obtiene mediante técnicas de elemento finito [Ci 1978] para
generación de las bases de los espacios solución de dimensión finita de los modelos.
Y finalmente para la implantación en computadora se requiere la aplicación de métodos numéricos
[Ci 1989] y lenguajes de programación.
37
2.1. Problema real.
Considere un cuerpo Ω homogéneo e isotrópico (en el sentido térmico), el cual es rígido y está fijo
en una región espacial, y el cual interactúa con su entorno de acuerdo las condiciones térmicas sobre
cada subfrontera Γi (i = 1, 2) , como se muestra en la figura 2.1, además considere que el cuerpo
alcanza un estado estacionario a partir de una condición inicial dada.
Fig. 2.1. Problema de Valores sobre la Frontera.
Revisando la figura observamos la existencia de tres regiones: el interior del cuerpo Ω y dos subfron-
teras Γi , i = 1, 2 que cumplen: Γ = ∪2i=1 Γi y además Ω = Ω ∪ Γ, donde Γ es la frontera del cuerpo
Ω.
Para cada región procedemos a analizar el problema, a través de las gráficas, correspondientes al
fenómeno de conducción de calor, entre el campo de la temperatura u y todos los valores posibles de
cantidad de calor q en cada punto x del cuerpo Ω.
En las subfronteras tenemos que el campo de flujo de calor q que pasa a través de la frontera puede
dirigirse hacia el interior o el exterior del cuerpo y su magnitud q (componente normal del campo de
flujo de calor q) se determina mediante −q = q · n, donde n es un campo vectorial normal unitario
exterior a la frontera Γ, notemos que el signo menos (−) corresponde a la convención: si −q es
positivo entonces q es saliente y cuando −q es negativo entonces q es entrante.
- Interior Ω.
En el interior del cuerpo el campo de temperatura (absoluta) u en cada punto x ∈ Ω está restringido
a valores positivos, y la cantidad de calor que hay en cada punto x ∈ Ω es a priori desconocida,
sin embargo consideremos una bola cerrada Br [x] de radio r > 0 centrada en x ∈ Ω a través de la
frontera de Br [x] , ∂Br [x] , entendemos que: hay un flujo de calor q (por conducción) que transporta
calor hacia el interior (+) o bien hacia el exterior (−) de Br [x] , existe además el calor generado o
consumido (fuente (+) o sumidero(−)) y el calor acumulado en la vecindad Br [x] deben cumplir con
la conservación de energía y por tanto se cumple con el balance de la cantidad de calor en Br [x] ,
38
gráficamente tenemos,
por lo que en este caso definimos −q como la suma (el balance) de la cantidad de calor (gráfica 1) y
éste cumple con el valor nulo.
Gráfica 1. Balance de Calor en el Interior.
- Subfrontera Γ1 .
En ésta subfrontera se tiene prescrita la temperatura y para cada punto x ∈ Γ1 el campo de temper-
atura u tiene como único valor admisible la temperatura prescrita û (x) , esto es, u (x) = û (x) , x ∈ Γ1
y para mantener la temperatura prescrita se tiene que la cantidad de calor −q (x) = q̂ (x) · n (x) que
pasa a través de cada punto x tiene como valores posibles: positivo, negativo o cero, esto es, para
mantener la temperatura prescrita es posible: inyectar calor (−q con valor negativo), extraer calor
(−q con valor positivo) o ni inyectar ni extraer (−q con valor nulo), y debido a que no podemos de-
terminar cual valor específico es requerido admitiremos como posibles todos los valores, en la gráfica
39
2 corresponden a la línea vertical.
Gráfica 2. Temperatura Prescrita.
- Subfrontera Γ2 .
En ésta subfrontera se tiene prescrita la cantidad de calor q̂ que pasa a través de cada punto x en
la frontera, donde −q̂ (x) = q̂ (x) · n (x) aquí q̂ es el flujo prescrito sobre Γ2 y para cada punto
x ∈ Γ2 el campo de temperatura u es restringido a valores positivos, esto es, la subfrontera Γ2 puede
alcanzar cualquier valor temperatura y para la temperatura alcanzada u (x) , x ∈ Γ2 se cumple que
−q (x) = −q̂ (x) , x ∈ Γ2 , esto es, sólo hay un valor posible para la cantidad de calor −q y coincide
con el valor prescrito −q̂, donde −q es positiva si se extrae calor y negativa si se inyecta calor, cuando
es nula la frontera está aislada térmicamente, en la gráfica 3 expresamos la condición.
Gráfica 3. Flujo Prescrito.
Con las gráficas anteriores en cada región se está tratando de interpretar que le sucede al cuerpo,
sin embargo debemos darnos cuenta de que estamos restringiendo o idealizando nuestro problema de
40
estudio, explícitamente se asume: el cuerpo es rígido (no hay cambios volumétricos por diferencia
de temperaturas), las fronteras son disyuntas y compatibles en sus condiciones sobre la frontera
(no se admite una zona de transición entre las condiciones), etc., e implícitamente se asume: no hay
cambios de fase del material, el material no se degrada en sus propiedades físicas, no hay condensación
o evaporación superficial, en fronteras de temperatura prescrita y temperatura mínima la cantidad
de calor inyectada o extraída es ilimitada, etc.
Por otra parte notemos que para la frontera de temperatura prescrita su gráfica correspondiente es
multivaluada.
Las gráficas anteriores es una forma de representar las ecuaciones de campo que forman el problema
de valores sobre la frontera (P V F ) , el cual es el modelo físico del problema de conducción de calor.
2.2. Modelo físico

Sea Ω ⊂ ξ ↔ R3 un cuerpo homogéneo e isotrópico (en el sentido térmico), el cual es rígido y está
fijo en una región espacial, y el cual interactua con su entorno de acuerdo las condiciones térmicas
sobre cada subfrontera Γi (i = 1, 2) , con frontera Γ = Γ1 ∪ Γ2 , donde Γ1 ∩ Γ2 = ∅, y considerando un
proceso de conducción de calor en el tiempo t ∈ [0, T ] ⊂ R expresamos el problema de valores sobre
la frontera y condiciones iniciales (P V F yCI), como:
Dada la condición inicial,
u (x, 0) = u0 en Ω. (2.1)
Encuentre u tal que cumpla:

d
−cv ρ u − div q + f = 0,
dt en Ω × [0, T ] , (2.2)
q = −K grad u,
u = û, sobre Γ1 × [0, T ] , (2.3)
−q · n = q̂, sobre Γ2 × [0, T ] , (2.4)
donde: q = −K grad u es la ley de Fourier de conducción de calor.
Identificamos los términos del modelo como:

u ≡ campo de temperatura [K] W ,
q ≡ campo de flujo de calor m 2 , W
div q ≡ divergencia del flujo de calor m 3 ,
f ≡ fuente de calor en el interior del cuerpo
W
W
m 3 ,
cv ρ dtd u ≡ término de acumulación m 3 ,
û ≡ temperatura prescrita [K] ,
q̂ ≡ cantidad de calor prescrita m W
2 ,

cv ≡ calor específico a volumen constante por unidad de masa KgJ K ,

ρ ≡ densidad de masa Kg m3
,
n ≡ campo de vectores normales unitarios exteriores a la frontera Γ,
41

K ≡ tensor de conductibilidad térmica mWK , K = kI en el caso isotrópico, aquí k es la constante
de conductibilidad térmica del material, I es el tensor identidad,
u0 ≡ temperatura inicial [K] en Ω para el tiempo t = 0.
Relacionando las ecuaciones (2.2 − 2.4) con las gráficas (1 − 4) respectivamente, se procede a rehacer
las gráficas:
Gráfica 4. Gráfica 5.
Gráfica 6.
Inicialmente procedemos a expresar el modelo físico [2.1 − 2.4] de acuerdo a ([Al 1989]) tenemos:
dada la condición inicial:
u (x, 0) = u0 (x) , x ∈ Ω, (2.5)

tal que se cumple:
d
−cv ρ u (x, t) − div (−K grad u (x, t)) + f (x, t) ∈ ∂φ (x, t; u (x, t)) = {0} , (2.6)
dt
si u (x, t) > 0, (x, t) ∈ Ω × [0, T ] .
42

 R, si γu (x, t) = û (x, t)
(−K grad u (x, t)) · n (x) ∈ ∂ψ 1 (x, t; γu (x, t)) = 0 < γu (x, t) < û (x, t) ,
 ∅, si (2.7)
û (x, t) < γu (x, t) ,
(x, t) ∈ Γ1 × [0, T ] .
(−K grad u (x, t)) · n (x) ∈ ∂ψ 2 (x, t; γu (x, t)) = {−q̂ (x, t)} ,
(2.8)
si γu (x, t) > 0, (x, t) ∈ Γ2 × [0, T ] .
Notar en gráficas (4 − 6) y en el modelo (2.5 − 2.8) se ha cambiado la temperatura u por la temprea-

tura γu donde γ representa el operador traza, el cual toma valores de la temperatura “muy cerca”
de la frontera y extiende los valores a la frontera. La extensión de valores hacia la frontera también
es necesaria cuando se requiere el cálculo de gradientes, los cuales se pueden obtener en el interior
del cuerpo pero no en la frontera.
2.3. Funciones primitivas.
Considere una función p de variable y valor real diferenciable. Grafiquemos la función p y su derivada
Dp donde se tiene que la derivada Dp tiene como valores la pendiente de la línea recta tangente a la
gráfica de la función p.
Gráfica de la función p. Gráfica de la función Dp.
43
Ahora considere la función Dp de variable y valor real y obtenga la primitiva correspondiente Dp dx,
esto es, la función Dp dx = q (x) + cte.
Gráfica de la función Dp. Gráfica de la función Dp dx.
Observemos que Dp dx define una familia de funciones q (x) + cte, x ∈ R, donde para cada valor
de la constante cte tenemos que q (x) + cte define una función afín, y cada funcion afín cumple
D (q + cte) = Dp.
En lo anterior se están manejando funciones de variable y valor real, donde se deben satisfacer
diversas condiciones para que los procesos de diferenciación e integración sean concretables, notando
además que las funciones utilizadas son univaluadas.
En las gráficas de nuestro problema se tiene funciones multivaluadas, y para su manejo vamos a
requerir ampliar conceptos como diferenciabilidad y trataremos que en el proceso de obtención de
primitivas, la selección de las constantes en el proceso de integración mantengan una interpretación
física para las primitivas.
Consideremos cada una de las gráficas 4 − 6 de nuestro modelo, y construyamos la gráfica de la
función primitiva correspondiente.
En el interior del cuerpo.
Gráfica 4.
44
La gráfica 4 representa lo que pasa en cada x ∈ Ω, la variable temperatura u está restringida a
u (x) > 0, los correspondientes valores representan el balance de la cantidad de calor en x ∈ Ω. Al
proceder a encontrar la primitiva correspondiente vemos que si en la gráfica 4 la variable temperatura
u (·) esdada en [K] y los valores del balance de calor −cv ρ dtd u (·) − div (−K grad u (·)) + f (·) es dado
en m W
3 tenemos que la primitiva ψ 1 es entonces dada en m W
3 K .
Gráfica 4a. Gráfica 4b.

Ahora cada valor del balance de la cantidad de calor (gráf ica 4) es la pendiente de líneas tangentes
a la primitiva, y son representadas en gráf ica 4a y que corresponden a una familia de funciones
afines. Y en la gráfica 4b se obtiene la gráfica de la función primitiva.
En la frontera de temperatura prescrita.
Gráfica 5.
La gráfica 5 representa lo que pasa en cada x ∈ Γ1 , la variable temperatura u está limitada a γu (x) =
u (x) , los correspondientes valores representan los posibles valores de cantidad de calor (aplicada)
sobre x ∈ Γ1 . Al proceder a encontrar la primitiva correspondiente vemos que si en la gráfica 5 la
45
variable temperatura
γu (·) es dado en [K] y los valores cantidad
W de calor (−K grad u (·)) · n (·) es
dado en m2 tenemos que la primitiva ψ 1 es entonces dada en m2 K .
W
Grafica 5a. Gráfica 5b.

Ahora, de los múltiples valores de cantidad de calor (gráf ica 5) cada uno es la pendiente de líneas
tangentes a la primitiva, y son representadas en gráf ica 5a y corta al eje horizontal en γθ (x) = θ (x) .
En la gráfica 5b se obtiene la gráfica de la función primitiva, la cual es la función que cae por arriba
o la misma altura de todas las líneas trazadas en la gráf ica 5a.
En la frontera de flujo prescrito.
Gráfica 6.
La gráfica 6 representa lo que pasa en cada x ∈ Γ2 , la variable temperatura u está restringida por
γu > 0 y los correspondientes valores representan los posibles valores de cantidad de calor (aplicada)
sobre x ∈ Γ2 . Al proceder a encontrar la primitiva correspondiente vemos que si en la gráfica 6 la
variable temperatura
γu (·) es dado en [K] y los valores cantidad
W decalor (−K grad θ (·)) · n (·) es
dado en m2 tenemos que la primitiva ψ 1 es entonces dada en m2 K .
W
46
Gráfica 6a. Gráfica 6b.
Ahora cada valor de la cantidad de de calor (gráf ica 6) cada uno es la pendiente de líneas tangentes
a la primitiva, y son representadas en gráf ica 6a, teniéndose una familia de funciones afines. y en
la gráfica 6b se obtiene la gráfica de la función primitiva buscada.
Y los correspondientes subpotenciales φ, ψ 1 y ψ 2 están dados por las expresiones siguientes:
φ (x, t; u (x, t)) = 0, u (x, t) > 0, (x, t) ∈ Ω × [0, T ] (2.9)


 0, γu (x, t) = û (x, t) > 0,
ψ 1 (x, t; γu (x, t)) = 0 < γu (x, t) < û (x, t) o , (x, t) ∈ Γ1 × [0, T ] . (2.10)
 +∞,
û (x, t) < γu (x, t) ,
ψ 2 (x, t; γu (x, t)) = − q̂ (x, t) γu (x, t) , γu (x, t) > 0, (x, t) ∈ Γ2 × [0, T ] . (2.11)
2.4. Modelos Matemáticos.

- Sea f : E −→ R ∪ {+∞} una función convexa (no necesariamente diferenciable), semi-continua
inferiormente s.c.i., con dominio efectivo, D (f ) ≡ {w ∈ E : f (w) < +∞} no vacío y cerrado. En-
tonces decimos que f es subdiferenciable en el punto u ∈ E, si existe al menos un u∗ ∈ E ′ tal
que,
f (w) ≥ f (u) + u∗ , w − uE ′ ×E ∀w ∈ E.
Esto es, f es subdiferenciable en u, si existe al menos una función afín w −→ f (u) + u∗ , w − uE ′ ×E ,
con pendiente u∗ que pasa através de f (u) , y cuya gráfica cae por debajo o a lo más a la misma
altura que f. El conjunto de los elementos u∗ ∈ E ′ (llamados subgradientes de f en u), define el
subdiferencial de f en u, denotado por ∂f (u) ; esto es,
∂f (u) = u∗ ∈ E ′ : f (w) ≥ f (u) + u∗ , w − uE ′ ×E ∀w ∈ E . (2.12)
47
La función primitiva de ∂f, f, es llamada el subpotencial de ∂f .
Observación. f : E −→ R ∪ {+∞} es convexa si y sólo si su subdiferencial ∂f : R −→ 2R es

monótono creciente.
Además, notemos que 2.12 puede ser expresada como:
u∗ ∈ ∂f (u) ⊂ E ′ ⇐⇒ f (w) ≥ f (u) + u∗ , w − uE′ ×E ∀w ∈ E
De lo expresado en la observación anterior, procedemos a modificar el modelo (2.5 − 2.8) y expresamos

las desigualdades variacionales locales correspondientes:
Dada la condición inicial:
u (x, 0) = u0 (x) , x ∈ Ω. (2.13)
Encuentre u tal que se cumpla:
u (x, t) ∈ D (φ (x, t; ·)) = (0, +∞) :
φ (x, t; η (x, t)) ≥ φ (x, t; u (x, t))
(2.14)
+ −cv ρ dtd u (x, t) − div (−K grad u (x, t)) + f (x, t) {η (x, t) − u (x, t)} ,
η (x, t) ∈ D (φ (x, t; ·)) , (x, t) ∈ Ω × [0, T ] .
γu (x, t) ∈ D (ψ 1 (x, t; ·)) = {û (x, t)} :

ψ 1 (x, t; γη (x, t)) ≥ ψ 1 (x, t; γu (x, t)) + {(−K grad u (x, t)) · n (x)} {γη (x, t) − γu (x, t)} , (2.15)
γη (x, t) ∈ D (ψ 1 (x, t; ·)) , (x, t) ∈ Γ1 × [0, T ] .
γu (x, t) ∈ D (ψ 2 (x; ·)) = (0, +∞) :
ψ 2 (x, t; γη (x, t)) ≥ ψ 2 (x, t; γu (x, t)) + {(−K grad u (x, t)) · n (x)} {γη (x, t) − γu (x, t)} , (2.16)
γη (x, t) ∈ D (ψ 2 (x, t; ·)) , (x, t) ∈ Γ2 × [0, T ] .
Sustituyendo los subpotenciales [2.9 − 2.11] y simplificando obtenemos:
u (x, 0) = u0 (x) , x ∈ Ω. (2.17)
Encuentre u tal que se cumpla:
u (x, t) ∈ D (φ (x, t; ·)) = (0, +∞) :
0 ≥ −cv ρ dtd u (x, t) − div (−K grad u (x, t)) + f (x, t) {η (x, t) − u (x, t)} , (2.18)
η (x, t) ∈ D (φ (x, t; ·)) , (x, t) ∈ Ω × [0, T ] .
γu (x, t) ∈ D (ψ 1 (x, t; ·)) = {û (x, t)} :

0 ≥ {(−K grad u (x, t)) · n (x)} {γη (x, t) − γu (x, t)} , (2.19)
γη (x, t) ∈ D (ψ 1 (x, t; ·)) , (x, t) ∈ Γ1 × [0, T ] .
γu (x, t) ∈ D (ψ 2 (x; ·)) = (0, +∞) :
{− q̂ (x, t)} {γη (x, t) − γu (x, t)} ≥ {(−K grad u (x, t)) · n (x)} {γη (x, t) − γu (x, t)} , (2.20)
γη (x, t) ∈ D (ψ 2 (x, t; ·)) , (x, t) ∈ Γ2 × [0, T ] .
48
Expresamos las desigualdades variacionales por región para el problema de conducción de calor.
u (x, 0) = u0 (x) , x ∈ Ω, (2.21)
Encuentre u tal que se cumple:
u (t)∈ D (Φ (t; ·)) :

0≥ −cv ρ dtd u (t) − div (−K grad u (t)) + f (t) {η (t) − u (t)} dΩ, (2.22)
Ω
η (t) ∈ D (Φ (t)) , t ∈ [0, T ] .
∈ D (Ψ1 (t; ·)) = {û (t)} :

γu (t)
0≥ {(−K grad u (t)) · n} {γη (t) − γu (t)} dΓ1 , (2.23)
Γ1
γη (t) ∈ D (Ψ1 (t)) , t ∈ [0, T ] .
γu (t) ∈ D (Ψ2 (t; ·)) :

{− q̂ (t)} {γη (t) − γu (t)} dΓ2 ≥ {(−K grad u (t)) · n} {γη (t) − γu (t)} dΓ2 , (2.24)
Γ2 Γ2
γη (x) ∈ D (Ψ2 (t)) , t ∈ [0, T ] .
los subpotenciales primales y dominios efectivos están dados por:


 Φ (t; η (t)) = φ (x, t; η (x, t)) dΩ = 0,
, (2.25)

 Ω
D (Φ (t; ·)) = {η : R → V (Ω) | η (t) > 0 en Ω} .


 Ψ1 (t; γη (t)) = ψ 1 (x, t; γη (x, t)) dΓ1 = 0,
, (2.26)

 Γ1
D (Ψ1 (t; ·)) = γη : R → B (Γ1 ) | γη (t)|Γ1 = û (t) sobre Γ1


 Ψ2 (t; γη (t)) = ψ 2 (x, t; γη (x, t)) dΓ2 = − q̂ (x, t) γη (x, t) dΓ2 ,
, (2.27)

 Γ2 Γ2
D (Ψ2 (t; ·)) = γη : R → B (Γ2 ) | γη (t)|Γ2 > 0 sobre Γ2
Ahora conjuntamos las desigualdades variacionales primales [2.22 − 2.24] , obteniéndose la formu-
lación variacional primal global fuerte.
u (x, 0) = u0 (x) , x ∈ Ω. (2.28)
49
Encuentre u (t) ∈ Kt :

d
cv ρ dt u (t) {η (t) − u (t)} dΩ + div (−K grad u (t)) {η (t) − u (t)} dΩ ≥
Ω Ω
f (t) {η (t) − u (t)} dΩ + q̂ (t) {γη (t) − γu (t)} dΓ2 + (2.29)
Ω Γ2
{(−K grad u (t)) · n} {γη (t) − γu (t)} dΓ, ∀η (t) ∈ Kt .

Γ=Γ1 ∪Γ2
Kt = D (Φ (t)) ∩ D
(Ψ1 (t) ◦ γ) ∩ D (Ψ2 (t) ◦ γ) (2.30)
= η : R → V Ω | η (t) > 0 en Ω, γη (t)|Γ1 = û (t) sobre Γ1 , γη (t)|Γ2 > 0 sobre Γ2 .
La formulación obtenida corresponde a una ecuación integro-diferencial, las derivadas parciales son
de segundo orden, para proceder a obtener una solución requerimos reducir el orden de las derivadas
parciales involucradas, para ello utilizaremos la fórmula de Green correspondiente al operador diver-
gencia.
Consideremos la relación div (α q) = α div q + q · grad α, integrando tenemos:

div (α q) dΩ = α div q dΩ + q · grad α dΩ
Ω Ω Ω
aplicando teorema de la divergencia al término izquierdo tenemos:

(α q) · n dΓ = α div q dΩ + q · grad α dΩ
Γ Ω Ω

(div q) α dΩ = − q · grad α dΩ + (q · n) α dΓ
Ω Ω Γ
sean q = −K grad u y α = {η − u} obtenemos la fórmula de Green correspondiente

div (−K grad u) {η − u} dΩ = (K grad u) · {grad η − grad u} dΩ (2.31)
Ω Ω

+ (−K grad u · n) {η − u} dΓ.
Γ
Aplicando la fórmula de Green [2.31] y simplificando obtenemos la formulación variacional global

primal débil.
u (0) = u0 en Ω, (2.32)
50
Encuentre u (t) ∈ Kt :

d
cv ρ dt u (t) {η (t) − u (t)} dΩ + K grad u (t) {grad η (t) − grad u (t)} dΩ ≥
Ω
Ω
(2.33)
f (t) {η (t) − u (t)} dΩ + q̂ (t) {γη (t) − γu (t)} dΓ2 , ∀η (t) ∈ Kt .
Ω Γ2

Kt = η : R → V Ω | η (t) > 0 en Ω, γη (t)|Γ1 = û (t) sobre Γ1 . (2.34)
2.5. Modelo Discreto.

El modelo matemático correspondiente a la formulación variacional débil, resulta ser el adecuado
para concretarse en un código computacional, por lo cual procedemos a expresar el modelo discreto
(espacial) así como la correspondiente versión coordenada.
Dada la condición inicial discreta:

uh (0) = u0 en Ωh . (2.35)
Encuentre uh (t) ∈ Kth :

cv ρ dtd uh (t) {η h (t) − uh (t)} dΩh + K grad uh (t) {grad ηh (t) − grad uh (t)} dΩh
Ωh Ωh
(2.36)
≥ f (t) {ηh (t) − uh (t)} dΩh + q̂ (t) {γη h (t) − γuh (t)} dΓ2h , ∀η h (t) ∈ Kth .
Ωh Γ2h

Kth = η h : R → Vh Ωh | η h (t) > 0 en Ωh , γη h (t)|Γ1h = û (t) sobre Γ1h . (2.37)
h
y donde se cumple que Vh Ωh , Kt es una familia que ”aproxima” a V Ω , Kt tal que:
1).- Vh Ωh es un espacio con producto interno de dimensión finita y

dim Vh Ωh = mh −→ +∞
h+↓0

2).- Kth es un subconjunto convexo de Vh Ωh no vacío,

3).- cv ρ dt uh (t) {ηh (t) − uh (t)} dΩh ”aproxima” a cv ρ dtd u (t) {η (t) − u (t)} dΩ.
d
Ω Ω
h
4).- K grad uh (t) {grad η h (t) − grad uh (t)} dΩh ”aproxima” a
Ωh

K grad u (t) {grad η (t) − grad u (t)} dΩ.
Ω
51

5).- f (t) {η h (t) − uh (t)} dΩh ”aproxima” a f (t) {η (t) − u (t)} dΩ
Ωh Ω

6).- q̂ (t) {γη h (t) − γuh (t)} dΓ2h ”aproxima” a q̂ (t) {γη (t) − γu (t)} dΓ2
Γ2h Γ2

7).- Vh Ωh ⊂ V Ω , aproximaciones internas.

Siendo Vh Ωh un espacio vectorial de dimensión finita, admitamos que una base de Vh Ωh sea

B = {w1 , w2 , · · · , wmh } , donde mh = dim Vh Ωh , entonces uh (t) ∈ Kth ⊂ Vh Ωh , admite la
representación:
uh (t) = υ j (t) wj (x) , j = 1, · · · , mh y η h (t) = hi (t) wi (x) , i = 1, · · · , mh .
Substituyendo y simplificando obtenemos la versión coordenada del modelo discreto.

υ (0) = υ 0 en Ωh , (2.38)
Encuentre υ ∈ Ktmh :
•
0 ≥ − M υ · {η − υ} − Aυ· {η − υ} + f · {η − υ} , ∀η ∈ Ktmh . (2.39)
Ktmh = η : R → Rmh | η (t) > 0 en Ωh , γη (t)|Γ1h = û sobre Γ1h . (2.40)

donde: 

 υ ∈ Ktmh , υ i = υ i (t) ,

 •

 υ j = dtd υ j (t) ,





 Mij = cv ρwj (x) wi (x) dΩh ,




Ω (2.41)

 Aij = K grad wj (x) grad wi (x) dΩh ,





 Ω





 fi = f (t) wi (x) dΩh + q̂ (t) γwi (x) dΓ2h .


Ω Γ2h

υ (0) = υ 0 en Ωh , (2.42)
Encuentre υ ∈ Ktmh :
•

IKm
t
h (η) ≥ IKmh (υ) +
t
− M υ −Aυ + f · {η − υ} , ∀η ∈ Ktmh . (2.43)
52
υ (0) = υ 0 en Ωh . (2.45)
•
Encuentre υ ∈ Ktmh : −M υ − Aυ + f ∈ ∂IKtmh (υ) . (2.46)
Nota. Aquí supondremos que toda solución de 2.46 tanto espacial como temporalmente siempre
satisface υ ∈ Ktmh con lo cual se cumplirá ∂IKm
t
h (υ) = {0} , ésta suposición para nuestros problemas
específicos se verifica numéricamente, encontrando los valores válidos de la fuente interna que hacen
que se satisfaga la suposición hecha.

υ (0) = υ 0 en Ωh . (2.48)
•
Encuentre υ ∈ Ktmh : −M υ − Aυ + f ∈ {0} . (2.49)
De acuerdo a la versión coordenada modelo semi-discreto espacial:

•
0 ∈ M υ +A (υ) − f. (2.51)
y donde expresamos la combinación convexa de A, con α ∈ (0, 1) ,

•
0 ∈ M υ +αA (υ) + (1 − α) A (υ) − f. (2.52)
donde:
υ ≡ vector coordenado de la temperatura [K],
•
υ≡ variación de la temperatura en el tiempo seg K
,
f ≡ el vector coordenado de fuente interna
W y acciones de frontera [W ] ,
A ≡ matriz de conductibilidad J térmica
K
,
M ≡ matriz de acumulación K ,
y para los algoritmos
W de integración temporal identificamos:
β = A K ,
J
µ = M K
,
ε ≡ error relativo normalizado K =[]
K
J
µ K
k ≡ paso de tiempo, k = f (ℓ) W = [seg] , donde f (ℓ) [ ] es una función que caracteriza a cada
β K
esquema de integración, de acuerdo al parámetro ℓ.
53
2.6. Método de Elemento Finito.
Consideremos el problema variacional lineal abstracto continuo,

encuentre u ∈ V Ω : (2.53)
a (u, v) = f (v) , ∀v ∈ V Ω .
El método de Galerkin, para la aproximación del problema (2.53)

consiste
en definir el problemas sim-
ilares en subespacios de dimensión
finita Vh Ωh del espacio V Ω , esto es, con cualquier subespacio
de dimensión finita Vh Ωh asociamos el problema

encuentre uh ∈ Vh Ωh : (2.54)
a (uh , vh ) = f (vh ) , ∀vh ∈ Vh Ωh .
Ahora,
para aplicar el método
de Garlekin, tenemos que construir el subespacio de dimensión finita
Vh Ωh del espacio V Ω .
Dado un modelo discreto para el cual su espacio solución Vh Ωh se puede representar por el espa-

cio generado por un conjunto de funciones (base), Vh Ωh = ϕ1 , ϕ2 , ϕ3 , . . . , ϕmh , con dimensión

dim Vh Ωh = mh , y la creación de una base particular permite concretar numéricamente las expre-
siones matriciales del modelo discreto y produce lo que llamaremos un Modelo Numérico.
En particular para éste propósito estudiaremos técnicas de elemento finito.

La idea a grandes rasgos es:
i).- Particionar el dominio espacial Ω en subdominios que se llamarán elementos finitos geométricos,
ii).- En cada elemento finito geométrico se definen los grados de libertad (los valores discretos) del
campo (variable) a conocer.
mh
iii).- Construir funciones base ϕj j=1 tales que ϕj valga 1 en el grado de libertad j − ésimo y 0 en
todos los demás.
Ver y estudiar en Ciarlet, P.G., The Finite Element Method for Elliptic Problems, North-Holland,
1978, el capítulo de elemento finitos (geométricos).
El método de elemento finito (mef) en su forma más simple, es un proceso específico para
construir los subespacios Vh Ωh , los cuales serán llamados espacios de elemento finito.
Esta construcción está caracterizada por tres aspectos básicos.
MEF 1.
Establecer una triangulación τ h sobre el dominio Ω, esto es, el cuerpo Ω es subdividido en un número
finito de subconjuntos Ej , llamados elemento finitos geométricos, tales que:
◦
a).- Ω ≈ Ωh ≡ Ej , E j = ∅.
Ej ∈τ h
◦ ◦
b).- E j E k = ∅, j = k.
54

c).- ∂Ej ∂Ek = ∅, si Ej y Ek son disjuntos ó ∂Ej ∂Ek es una cara, arista o vértice común si Ej
y Ek son adyacentes.
Una vez establecida la triangulación τ h sobre Ω, se define
Xh ⊂ f : Ω → R ,
(2.55)
Xh = [wj ]m
j=1 , dim Xh = mh ,
h
a este espacio se le conoce como Espacio de Elemento Finito Global .

Dado un Espacio de Elemento Finito Global se definen los espacios
PE = { vh |E : vh ∈ Xh , E ∈ τ h } (2.56)
estos espacios son llamados Espacios

de Elemento Finito Local .
Consecuentemente el espacio Vh Ωh viene a ser

Vh Ωh = {vh ∈ Xh : vh satisfaga condiciones de frontera} (2.57)
M E F 2.
Los espacios PE , E ∈ τ h , son espacios polinomiales o funciones cercanas a polinomios .
M E F 3.
Existe una base canónica en el espacio discreto Vh Ωh cuyas funciones base correspondiente tiene
soportes compactos, los cuales son tan pequeño como sea posible, siendo implícitamente entendido
que estas funciones base son fácilmente descritas.
-Por Elemento Finito entendemos a la terna
(E, PE , ΣE ) (2.58)
donde
E ≡ el elemento finito geométrico,
PE ≡ el espacio de elemento finito local, y
ΣE ≡ el conjunto de grados de libertad
siendo
ΣE = {φi ∈ L (PE , R) : φi sea PE -unisolvente} (2.59)
además
0, si i = j
φi (p (xj )) = δ ij = , p (·) ∈ PE , xj ∈ E. (2.60)
1, si i = j
55
2.6.1. Elementos unidimensionales.
Caso lineal escalar.

Sea E el elemento finito geométrico
el espacio de elemento finito local PE es
PE = {p : E → R, p (x) = β 0 + β 1 x}
y el conjunto de grados de libertad ΣE es
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE }
Construyamos el espacio de elemento finito local, PE , esto es, construimos las funciones base que
generen dicho espacio, PE = [w̃1 , w̃2 ] ,
PE = {v : E → R | v (x) = α1 w̃1 (x) + α2 w̃2 (x)} ,
tenemos de la propiedad de PE -unisolvencia que:
φi (p (xj )) = δ ij entonces w̃i (x) = φi (p) (x) , i = 1, 2.
para la primera función base tenemos:

φ1 (p (b1 )) = β 0 + β 1 b1 = 1 1 b1 β0 1
⇒ = ⇒
φ1 (p

(b2 )) =
β 0 + β 1 b2 = 0 1 b2 β1 0

1 b1 1 1

0 b2 b2 1 0 −1 b2 − x
β 0 = =
y β 1 =

=
=⇒ w̃1 (x) = , x ∈ E,
1 b1 b2 − b1 1 b1 b2 − b1 b2 − b1

1 b2 1 b2
para la segunda función base tenemos:
φ2 (p (b1 )) = β 0 + β 1 b1 = 0 1 b1 β0 0
⇒ =
φ2 (p (b 2 )) = β 0 + β 1 b2 = 1
1 b2 β1 1

0 b1 1 0

1 b2 −b1 1 1 1 x − b1
⇒ β0 =
=
y β1 = =
=⇒ w̃2 (x) = , x ∈ E,
1 b1 b2 − b1 1 b1 b2 − b1 b2 − b1

1 b2 1 b2
56
gráficamente tenemos:
Nota. Los coeficientes de las dos funciones basew̃1 (x) , w̃2 (x) pueden ser obtenidas en un solo
1 b1
proceso calculando la inversa de la matriz , esto es:
1 b2
−1
1 b1 − b1b−b
2 b1
= 1
2 b1 −b
1
2 , donde los coeficentes están dados por columnas.
1 b2 b1 −b2
− b1 −b 2
Caso cuadrático escalar.

PE = p : E → R, p (x) = β 0 + β 1 x + β 2 x2
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE , φ3 (p) ∈ PE } .
generen dicho espacio, PE = [w̃1 , w̃2 , w̃3 ] ,
PE = {v : E → R | v (x) = α1 w̃1 (x) + α2 w̃2 (x) + α3 w̃3 (x)} ,
φi (p (xj )) = δ ij entonces w̃i (x) = φi (p) (x) , i = 1, · · · , 3.
para la primera función base tenemos:
57
 2     
 φ1 (p (b1 )) = β 0 + β 1 b1 + β 2 (b1 ) = 1 1 b1 (b1 )2  β 0   1 
φ1 (p (b2 )) = β 0 + β 1 b2 + β 2 (b2 )2 = 0 ⇒  1 b2 (b2 )2  β 1 = 0 ⇒
 2 2    
φ1 (p (b3 )) = β 0 + β 1 b3 + β 2 (b3 ) = 0 1 b3 (b3 ) β 2 0
2

1 b1 (b1 )

2
0 b2 (b2 )

0 b3 (b3 )2 b2 b3 (b3 − b2 ) b2 b3
β 0 = 2 = 2 2 2 2 2 2
= 2
,
1 b1 (b1 )

b2 b 3 − b 2 b3 − b 1 b 3 + b 1 b3 + b1 b2 − b 1 b 2 b 2 b 3 − b2 b1 + b1 − b 3 b 1
2
1 b2 (b2 )

1 b3 (b3 )2
2

1 1 (b1 )

2
1 0 (b2 )

1 0 (b3 )2 b22 − b23 − (b2 + b3 )
β 1 = 2 = 2 2 2 2 2 2
= ,y
1 b1 (b1 )

b2 b3 − b2 b3 − b1 b3 + b1 b3 + b1 b2 − b1 b2 b2 b3 − b2 b1 + b21 − b3 b1
2
1 b2 (b2 )

1 b3 (b3 )2

1 b1 1

1 b2 0

1 b3 0 b3 − b2 1
β 2 = 2 = 2 2 2 2 2 2
= ,
1 b1 (b1 )

b2 b3 − b2 b3 − b1 b3 + b1 b3 + b1 b2 − b1 b2 b2 b3 − b2 b1 + b21 − b3 b1
2
1 b2 (b2 )

1 b3 (b3 )2
b2 b3 − (b2 + b3 ) x + x2
w̃1 (x) = , x ∈ E,
b2 b3 − (b2 + b3 ) b1 + b21
para
 la segunda función base tenemos:     
2
 φ2 (p (b1 )) = β 0 + β 1 b1 + β 2 (b1 ) = 0 1 b1 (b1 )2  β 0   0 
φ (p (b2 )) = β 0 + β 1 b2 + β 2 (b2 )2 = 1 ⇒  1 b2 (b2 )2  β 1 = 1 ⇒
 2 2 2    
φ2 (p (b3 )) = β 0 + β 1 b3 + β 2 (b3 ) = 0 1 b3 (b3 ) β2 0
2

0 b1 (b1 )

2
1 b2 (b2 )

0 b3 (b3 )2 − (b3 − b1 ) b1 b3 −b1 b3
β 0 = 2 = 2 2 2 2 2 2
= ,
1 b1 (b1 )

b2 b3 − b2 b3 − b1 b3 + b1 b3 + b1 b2 − b1 b2 −b3 b1 + b2 b1 − b22 + b2 b3
2
1 b2 (b2 )

1 b3 (b3 )2
2

1 0 (b1 )

2
1 1 (b2 )

1 0 (b3 )2 b23 − b21 (b1 + b3 )
β 1 = 2 = 2 2 2 2 2 2
= 2
,y
1 b1 (b1 )

b2 b 3 − b 2 b3 − b 1 b 3 + b 1 b3 + b1 b2 − b 1 b 2 −b3 b1 + b2 b1 − b 2 + b 2 b 3
2
1 b2 (b2 )

1 b3 (b3 )2
58

1 b1 0

1 b2 1

1 b3 0 − (b3 − b1 ) −1
β2 = = = ,

1 b1 (b1 )2

b2 b23 − b22 b3 2 2 2 2
− b1 b3 + b1 b3 + b1 b2 − b1 b2 −b3 b1 + b2 b1 − b22 + b2 b3

1 b2 (b2 )2

1 b3 (b3 )2
−b3 b1 + (b1 + b3 ) x − x2
w̃2 (x) = , x ∈ E,
−b3 b1 + (b1 + b3 ) b2 − b22
para
 la tercera función base tenemos:2     
 φ3 (p (b1 )) = β 0 + β 1 b1 + β 2 (b1 ) = 0 1 b1 (b1 )2  β 0   0 
φ (p (b2 )) = β 0 + β 1 b2 + β 2 (b2 )2 = 0 ⇒  1 b2 (b2 )2  β 1 = 0 ⇒
 3 2 2    
φ3 (p (b3 )) = β 0 + β 1 b3 + β 2 (b3 ) = 1 1 b3 (b3 ) β2 1
2

0 b1 (b1 )

2
0 b2 (b2 )

1 b3 (b3 )2 b1 b2 (b2 − b1 ) b2 b3
β 0 = 2 = 2 2 2 2 2 2
= ,
1 b1 (b1 )

b2 b3 − b2 b3 − b1 b3 + b1 b3 + b1 b2 − b1 b2 −b3 b1 + b2 b1 + b23 − b2 b3
2
1 b2 (b2 )

1 b3 (b3 )2
2

1 0 (b1 )

2
1 0 (b2 )

1 1 (b3 )2 b21 − b22 − (b1 + b2 )
β 1 = 2 = 2 2 2 2 2 2
= 2
,y
1 b1 (b1 )

b2 b 3 − b 2 b3 − b 1 b 3 + b 1 b3 + b1 b2 − b 1 b 2 −b3 b1 + b2 b1 + b3 − b 2 b 3
2
1 b2 (b2 )

1 b3 (b3 )2

1 b1 0

1 b2 0

1 b3 1 b2 − b1 1
β 2 = 2 = 2 2 2 2 2 2
= 2
,
1 b1 (b1 )

b2 b 3 − b 2 b3 − b 1 b 3 + b 1 b3 + b1 b2 − b 1 b 2 −b3 b1 + b2 b1 + b3 − b 2 b 3
2
1 b2 (b2 )

1 b3 (b3 )2
b2 b1 − (b1 + b2 ) x + x2
w̃3 (x) = , x ∈ E,
b2 b1 − (b1 + b2 ) b3 + b23
59
Nota. Los coeficientes de las tres funciones base  w̃1 (x) , w̃22(x)  , w̃3 (x) pueden ser obtenidas en un
1 b1 (b1 )
solo proceso calculando la inversa de la matriz 1 b2 (b2 )2  , esto es:

1 b3 (b3 )2
  
2 −1 b2 b3 −b3 b1 b2 b1
1 b1 (b1 ) 2
b2 b3 −b2 b1 +b1 −b3 b1 2
−b3 b1 +b2 b1 −b2 +b2 b3 2
−b3 b1 +b2 b1 +b3 −b2 b3
 1 b2 (b2 )2  =  − b2 +b3 b1 +b3
 b2 b3 −b2 b1 +b1 −b3 b1 −b3 b1 +b2 b1 −b2 +b2 b3
2 2 − b1 +b2 
−b3 b1 +b2 b1 +b23 −b2 b3  , donde los coe-
1 b3 (b3 )2 1
b2 b3 −b2 b1 +b21 −b3 b1
− −b3 b1 +b2 b11 −b2 +b2 b3 −b3 b1 +b2 b11 +b2 −b2 b3
2 3
ficientes están dados por columnas.
Caso cúbico escalar.

PE = p : E → R, p (x) = β 0 + β 1 x + β 2 x2 + β 3 x3
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE , φ3 (p) ∈ PE , φ4 (p) ∈ PE } .
generen dicho espacio, PE = [w̃1 , w̃2 , w̃3 , w̃4 ] ,
PE = {v : E → R | v (x) = α1 w̃1 (x) + α2 w̃2 (x) + α3 w̃3 (x) + α4 w̃4 (x)} ,
φi (p (xj )) = δ ij entonces w̃i (x) = φi (p) (x) , i = 1, · · · , 4.
60
para este elemento finito geométrico procederemos a calcular los coeficientes de las cuatro funciones
base
 en un solo proceso:

 φ1−4 (p (b1 )) = β 0 + β 1 b1 + β 2 (b1 )2 + β 3 (b1 )3 = 1, 0, 0, 0

φ1−4 (p (b2 )) = β 0 + β 1 b2 + β 2 (b2 )2 + β 3 (b2 )3 = 0, 1, 0, 0
⇒

 φ1−4 (p (b3 )) = β 0 + β 1 b3 + β 2 (b3 )2 + β 3 (b3 )3 = 0, 0, 1, 0

φ1−4 (p (b4 )) = β 0 + β 1 b4 + β 2 (b4 )2 + β 3 (b4 )3 = 0, 0, 0, 1
          
1 b1 (b1 )2 (b1 )3   β0     1 
  0 
   0     0  
 1 b2 (b2 )2 (b2 )3   β 1   0   1   0   0 
 
 1 b3 (b3 )2 (b3 )3   β 2  =  0  ,  0  ,  1  ,  0  ⇒

 
               
1 b4 (b4 )2 (b4 )3 β3 0 0 0 1
   2 3 −1    −b2 b3 b4 
 β  1 b 1 (b 1 ) (b 1 )  1    b b b −b b b −b2 b +b b b −b2 b +b b b +b3 −b2 b 


 0    2 3 

    
2 3 1 2 3 4 1 2 2 4 1 1 3
b3 b2 +b2 b4 +b3 b4
3 4 1 1 1 4 

β1  1 b 2 (b 2 ) (b 2 )  0 2 2 3 2
b2 b3 b1 −b2 b3 b4 −b1 b2 +b2 b4 b1 −b1 b3 +b3 b4 b1 +b1 −b1 b4
= =

 β  1 b3 (b3 )2 (b3 )3   0   b2 +b3 +b4
− b2 b3 b1 −b2 b3 b4 −b2 b2 +b 
 2   2 3

    
 1 2 b4 b1 −b1 b3 +b3 b4 b1 +b1 −b1 b4 
2 3 2


β3 1 b4 (b4 ) (b4 ) 0 1
b b b −b b b −b2 b +b b b −b2 b +b b b +b3 −b2 b
2 3 1 2 3 4 1 2 2 4 1 1 3 3 4 1 1 1 4
−b2 b3 b4 + (b3 b2 + b2 b4 + b3 b4 ) x − (b2 + b3 + b4 ) x2 + x3

w̃1 (x) = , x ∈ E,
−b2 b3 b4 + (b2 b3 + b2 b4 + b3 b4 ) b1 − (b2 + b3 + b4 ) b21 + b31
   −1    b3 b4 b1 

 β0 
 1 b1 (b1 )2 (b1 )3 
 0 


 −b2 b3 b1 +b3 b4 b1 +b1 b22 −b2 b4 b1 +b22 b3 −b2 b3 b4 −b32 +b22 b4 


      − b3 b1 +b4 b1 +b3 b4 
β1 1 b2 (b2 )2 (b2 )3  1 −b2 b3 b1 +b3 b4 b1 +b1 b22 −b2 b4 b1 +b22 b3 −b2 b3 b4 −b32 +b22 b4
=  =
  
 β2  1 b3 (b3 )2 (b3 )3    0  

b1 +b3 +b4
−b2 b3 b1 +b3 b4 b1 +b1 b22 −b2 b4 b1 +b22 b3 −b2 b3 b4 −b32 +b22 b4


    
 − 

β3 1 b4 (b4 )2 (b4 )3 0 1
−b2 b3 b1 +b3 b4 b1 +b1 b2 −b2 b4 b1 +b2 b3 −b2 b3 b4 −b3 +b2 b4
2 2 2 2
b3 b4 b1 − (b3 b1 + b4 b1 + b3 b4 ) x + (b1 + b3 + b4 ) x2 − x3
w̃2 (x) = , x ∈ E,
b3 b4 b1 − (b3 b1 + b4 b1 + b3 b4 ) b2 + (b1 + b3 + b4 ) b22 − b32
   −1    −b2 b4 b1 

 β 
 1 b1 (b1 )2 (b1 )3 
 0  
 −b2 b4 b1 +b2 b3 b1 +b3 b4 b1 −b1 b23 +b2 b3 b4 −b2 b23 −b23 b4 +b33 

 0
       b1 b2 +b4 b1 +b2 b4 

β1 1 b2 (b2 )2 (b2 )3  0 −b2 b4 b1 +b2 b3 b1 +b3 b4 b1 −b1 b23 +b2 b3 b4 −b2 b23 −b23 b4 +b33
=

 =

 β  1 b3 (b3 )2 (b3 )3   1   − −b2 b4 b1 +b2 b3 b1 +b3 b4 bb11 +b 2 +b4

 2   2 3

    

−b1 b23 +b2 b3 b4 −b2 b23 −b23 b4 +b33 


β3 1 b4 (b4 ) (b4 ) 0 1
−b b b +b b b +b b b −b b2 +b b b −b b2 −b2 b +b3
2 4 1 2 3 1 3 4 1 1 3 2 3 4 2 3 3 4 3
−b2 b4 b1 + (b1 b2 + b4 b1 + b2 b4 ) x − (b1 + b2 + b4 ) x2 + x3

w̃3 (x) = , x ∈ E,
−b2 b4 b1 + (b2 b1 + b4 b1 + b2 b4 ) b3 − (b1 + b2 + b4 ) b23 + b33
   −1    b2 b3 b1 

 β0 
 1 b1 (b1 )2 (b1 )3 
 0  −b2 b4 b1 +b2 b3 b1 +b1 b24 −b3 b4 b1 +b2 b24 −b2 b3 b4 −b34 +b3 b24
      − b1 b2 +b3 b1 +b3 b2 
β1 1 b2 (b2 )2 (b2 )3  0  −b2 b4 b1 +b2 b3 b1 +b1 b24 −b3 b4 b1 +b2 b24 −b2 b3 b4 −b34 +b3 b24 
=  = 
  
 β2  1 b3 (b3 )2 (b3 )3    0   −b2 b4 b1 +b2 b3 b1 +b1 b2 b−b 1 +b2 +b3
2 3
3 b4 b1 +b2 b4 −b2 b3 b4 −b4 +b3 b4
2 
    4
β3 1 b4 (b4 )2 (b4 )3 1 1
− −b2 b4 b1 +b2 b3 b1 +b1 b2 −b3 b4 b1 +b2 b2 −b2 b3 b4 −b3 +b3 b2
4 4 4 4
b2 b3 b1 − (b1 b2 + b3 b1 + b3 b2 ) x + (b1 + b2 + b3 ) x2 − x3
w̃4 (x) = , x ∈ E,
b2 b3 b1 − (b2 b1 + b3 b1 + b3 b2 ) b4 + (b1 + b2 + b3 ) b24 − b34
61
2.6.2. Elementos bidimensionales.

Caso triángulo escalar.
PE = {p : E → R, p (x, y) = β 0 + β 1 x + β 2 y}
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE , φ3 (p) ∈ PE } .
generen dicho espacio, PE = [w̃1 , w̃2 , w̃3 ] ,
PE = {v : E → R | v (x, y) = α1 w̃1 (x, y) + α2 w̃2 (x, y) + α3 w̃3 (x, y)} ,
φi (p (xj , yj )) = δ ij entonces w̃i (x, y) = φi (p) (x, y) , i = 1, · · · , 3.
para este elemento finito geométrico procederemos a calcular los coeficientes de las tres funciones
base
 φ1−3 (p (bx1 , by1 )) = β 0 + β 1 bx1 + β 2 by1 = 1, 0, 0
φ (p (bx2 , by2 )) = β 0 + β 1 bx2 + β 2 by2 = 0, 1, 0 ⇒
 1−3
 φ1−3 (p (bx3 ,by
3 )) = β 0 +β 1 b
x3 +
β 2 b
y3 =
0, 0, 1
1 bx1 by1  β 0   1   0   0 
 1 bx2 by2  β 1 = 0 , 1 , 0 ⇒
       
1 bx3 by3 β2 0 0 1
62
 
   −1    bx2 by3 −by2 bx3

 β0  1 bx1 by1  1   bx2 by3 −by2 bx3 −bx1 by3 +by1 bx3 +bx1 by2 −by1 bx2 

−by3 +by2
β1 =  1 bx2 by2  0 = bx2 by3 −by2 bx3 −bx1 by3 +by1 bx3 +bx1 by2 −by1 bx2
      

β2 1 bx3 by3 0  − −bx3 +bx2 
bx by −by bx −bx by +by bx +bx by −by bx
2 3 2 3 1 3 1 3 1 2 1 2
(bx2 by3 − by2 bx3 ) + (−by3 + by2 ) x − (−bx3 + bx2 ) y

w̃1 (x, y) = , (x, y) ∈ E,
(bx2 by3 − by2 bx3 ) + (−by3 + by2 ) bx1 − (−bx3 + bx2 ) by1
 
   −1    − bx1 by3 −by1 bx3

 β0  1 bx1 by1  0   bx2 by3 −by2 bx3 −bx1 by3 +by1 bx3 +bx1 by2 −by1 bx2 
 
−by3 +by1
β 
= 1 bx2 by2  1 = −
 1     
bx2 by3 −by2 bx3 −bx1 by3 +by1 bx3 +bx1 by2 −by1 bx2


β2 1 bx3 by3 0  −bx3 +bx1 
− (bx1 by3 − by1 bx3 ) − (−by3 + by1 ) x + (−bx3 + bx1 ) y

w̃2 (x, y) = , (x, y) ∈ E,
− (bx1 by3 − by1 bx3 ) − (−by3 + by1 ) bx2 + (−bx3 + bx1 ) by2
 
   −1    bx1 by2 −by1 bx2

 β0  1 bx1 by1  0    bx2 by3 −by2 bx3 −bx1 by3 +by1 bx3 +bx1 by2 −by1 bx2  
−by2 +by1
β 
= 1 bx2 by2  0 =
 1     


β2 1 bx3 by3 1  − −bx2 +bx1 
(bx1 by2 − by1 bx2 ) + (−by2 + by1 ) x − (−bx2 + bx1 ) y

w̃3 (x, y) = , (x, y) ∈ E,
(bx1 by2 − by1 bx2 ) + (−by2 + by1 ) bx3 − (−bx2 + bx1 ) by3
Caso triángulo escalar 2.

63
PE = p : E → R, p (x, y) = β 0 + β 1 x + β 2 y + β 3 x2 + β 4 xy + β 5 y 2
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE , φ3 (p) ∈ PE , φ4 (p) ∈ PE , φ5 (p) ∈ PE , φ6 (p) ∈ PE } .
TAREA 1. Construya el espacio de elemento finito local, PE , esto es, construimos las funciones
base que generen dicho espacio, PE = [w̃1 , w̃2 , w̃3 , w̃4 , w̃5 , w̃6 ] .

Caso cuadrángulo escalar.
PE = {p : E → R, p (x, y) = β 0 + β 1 x + β 2 y + β 3 xy}
generen dicho espacio, PE = [w̃1 , w̃2 , w̃3 , w̃4 ] ,
PE = {v : E → R | v (x, y) = α1 w̃1 (x, y) + α2 w̃2 (x, y) + α3 w̃3 (x, y) + α4 w̃4 (x, y)} ,
φi (p (xj , yj )) = δ ij entonces w̃i (x, y) = φi (p) (x, y) , i = 1, · · · , 4.
para este elemento finito geométrico procederemos a calcular los coeficientes de las cuatro funciones
base

 φ (p (bx1 , by1 )) = β 0 + β 1 bx1 + β 2 by1 + β 3 bx1 by1 = 1, 0, 0, 0
 1−4
φ1−4 (p (bx2 , by2 )) = β 0 + β 1 bx2 + β 2 by2 + β 3 bx2 by2 = 0, 1, 0, 0
⇒

 φ1−4 (p (bx3 , by3 )) = β 0 + β 1 bx3 + β 2 by3 + β 3 bx3 by3 = 0, 0, 1, 0

φ (p (bx4 , by4 )) = β0 + β 1bx4 + β 2 by4 +
β 3 bx4 b
y4 = 
0, 0,0, 1 
 1−4 
1 bx1 by1 bx1 by1   β0     1     0     0    0  
 1 bx2 by2 bx2 by2   β 1   0   1   0   0 
 
 1 bx3 by3 bx3 by3   β 2  =  0  ,  0  ,  1  ,  0 

 
                
1 bx4 by4 bx4 by4 β3 0 0 0 1
64
Haciendo
 las consideraciones
  geométricas
  bx3= bx 1
, b
x4 = b by2 =by1 , by4 = by3 entonces
x2 , 
1 bx1 by1 bx1 by1   β0     1 
  0 
   0     0  
 1 bx2 by1 bx2 by1   β 1   0   1   0   0 
 
 1 bx1 by3 bx1 by3   β 2  =  0  ,  0  ,  1  ,  0  ⇒

 
               
1 bx2 by3 bx2 by3 β3 0 0 0 1
 
   −1    bx2 by3

 β0  1 bx1 by1 bx1 by1  1  
 −bx1 by3 +bx2 by3 −bx2 by1 +bx1 by1  

 
   
   −  b 

β1 1 b b b b  0
y3
=  x2 y1 x2 y1 
= −bx1 by3 +bx2 by3 −bx2 by1 +bx1 by1

 β  1 bx1 by3 bx1 by3   0   − −bx by +bx by x−b
b 2

 2    
   


x2 by1 +bx1 by1 

β3 1 bx2 by3 bx2 by3 0  1 3 2
1
3

−bx1 by3 +bx2 by3 −bx2 by1 +bx1 by1
(bx2 by3 ) − (by3 ) x − (bx2 ) y + xy

w̃1 (x, y) = , (x, y) ∈ E,
bx2 by3 − by3 bx1 − bx2 by1 + bx1 by1
 
   −1    − bx1 by3

β 1 bx1 by1 b x1 by1  
   0 
   
−b x by +b x2 by3 −bx2 by1 +bx1 by1

 0 
   1 3

β1 1 bx2 by1 b x2 by1   1   −b b +b b by−b 3 
=


  0   = x1 y3 x2 y3
bx1
x2 by1 +bx1 by1

 β 
 1 bx1 by3 b x1 by3  
 2
     
 −bx by +bx by −bx by +bx by 


β3 1 bx2 by3 b x2 by3 0  − 1 3 2 31 2 1 1 1 
−bx by +bx by −bx by +bx by
1 3 2 3 2 1 1 1
− (bx1 by3 ) + (by3 ) x + (bx1 ) y − xy

w̃2 (x, y) = , (x, y) ∈ E,
−bx1 by3 + by3 bx2 + bx1 by1 − bx2 by1
 
   −1    − bx2 by1

 β  1 bx1 by1 b x1 by1  0  
    
−b x by +b x2 by3 −bx2 by1 +bx1 by1

 0 
   1 3

β1 1 bx2 by1 b x2 by1   0   −b b +b b by−b 1 
=


  1  = x1 y3 x2 y3
bx2
x2 by1 +bx1 by1

 β 
 1 bx1 by3 b x1 by3  
 2
    
 −bx by +bx by −bx by +bx by 


β3 1 bx2 by3 b x2 by3 0  − 1 3 2 31 2 1 1 1 
−bx by +bx by −bx by +bx by
1 3 2 3 2 1 1 1
− (bx2 by1 ) + (by1 ) x + (bx2 ) y − xy

w̃3 (x, y) = , (x, y) ∈ E,
 
   −1    bx1 by1

 β  1 bx1 by1 b x1 by1  0  
    
−bx b y +bx 2 by3 −bx2 by1 +bx1 by1

 0 
   1 3

β1 1 bx2 by1 b x2 by1   0   − −b b +b b by−b 1 
=

 =
  0   −
x1 y3 x2 y3
bx1
x2 by1 +bx1 by1

 β 
 1 bx1 by3 b x1 by3  
 2
    
 −bx1 by3 +bx2 by3 −bx2 by1 +bx1 by1 


β3 1 bx2 by3 b x2 by3 1  1 
−bx1 by3 +bx2 by3 −bx2 by1 +bx1 by1
(bx1 by1 ) − (by1 ) x − (bx1 ) y + xy

w̃4 (x, y) = , (x, y) ∈ E,
65
Caso cuadrángulo escalar 2.


PE = p : E → R, p (x, y) = β 0 + β 1 x + β 2 y + β 3 x2 + β 4 xy + β 5 y 2 + β 6 x2 y + β 7 xy 2 + β 8 x2 y 2
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE , · · · , φ9 (p) ∈ PE } .
2.6.3. Elementos tridimensionales.

Caso tetrahedro.

PE = {p : E → R, p (x, y, z) = β 0 + β 1 x + β 2 y + β 3 z}
66
TAREA 2. Construya el espacio de elemento finito local, PE , esto es, construimos las funciones base
que generen dicho espacio, PE = [w̃1 , w̃2 , w̃3 , w̃4 ] .

Caso pentahedro.
PE = {p : E → R, p (x, y, z) = β 0 + β 1 x + β 2 y + β 3 z + β 4 yz + β 5 xz}
ΣE = {φi ∈ L (PE , R) : φ1 (p) , φ2 (p) , φ3 (p) , φ4 (p) , φ5 (p) , φ6 (p) ∈ PE } .
Caso hexahedro.
PE = {p : E → R, p (x, y, z) = β 0 + β 1 x + β 2 y + β 3 z + β 4 xy + β 5 yz + β 6 xz + β 7 xyz}
ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE , · · · , φ8 (p) ∈ PE } .
67
2.7. Conducción de Calor Unidimensional.
Regresando al problema de conducción de calor evolutivo, consideremos el cuerpo Ω con las dimen-
siones geométricas:
Fig. 2.2. Geometría problema ejemplo.
bajo una discretización espacial:
Fig. 2.3. Triangulación.
consideremos el elemento finito local:
Fig. 2.4. Elemento Finito Geométrico.
con los valores de las propiedades físicas y dimensiones siguientes:

cv ≡ calor específico a volumen constante por unidad de masa KgJ K , cv = 480.

, ρ = 7850.
A ≡ área transversal [m2 ] , A = 0.001.
long = L ≡ longitud del cuerpo [m] , long = 1.0.
n ≡ número de elementos finitos de la discretización.
long
b2 − b1 ≡ longitud de los elementos finitos [m] , b2 − b1 = .
W n.
k ≡ constante de conductibilidad
W térmica m K , k = 50.
f ≡ fuente interna m 3 , f = 10000.
De acuerdo a 2.41 tenemos la definición de los componentes de las matrices y vectores del modelo en
estudio, para el elemento finito local explicitamos las versiones locales correspondientes obteniéndose:
68
L L
Mij = cv ρwj (x) wi (x) dΩh = cv ρwj (x) wi (x) dA dx = cv ρAwj (x) wi (x) dx
Ωh 0 A 0
a2 a
n+1
= cv ρAwj (x) wi (x) dx + · · · + cv ρAwj (x) wi (x) dx

a1 an

= cv ρA wj (x)|E1 wi (x)|E1 dx + · · · + cv ρA wj (x)|En wi (x)|En dx
E1 En  
#n n
# n
#
= cv ρA wj (x)|Ek wi (x)|Ek dx =  cv ρAŵm (x) ŵl (x) dx = Ek
Mlm
(j=m,i=l)
k=1 E k=1 Ek k=1
k (j=m,i=l)

donde MEk
lm = cv ρAŵm (x) ŵl (x) dx es la matriz de acumulación local y la expresión
Ek
n
#
MEk
lm corresponde al proceso conocido como ensamble.
(j=m,i=l)
k=1
L
Aij = K grad wj (x) grad wi (x) dΩh = K grad wj (x) grad wi (x) dA dx
Ωh 0 A
L
= AK grad wj (x) grad wi (x) dx
0
a2 a
n+1
= AK grad wj (x) grad wi (x) dx + · · · + AK grad wj (x) grad wi (x) dx

a1 an

= AK grad wj (x)|E1 grad wi (x)|E1 dx + · · · + AKgrad wj (x)|En grad wi (x)|En dx
E1 En
n
#
= AKgrad wj (x)|Ek grad wi (x)|Ek dx
k=1 E
k 
n
# n
#
=  AKgradŵm (x) gradŵl (x) dx = Ek
Alm
(j=m,i=l)
k=1 Ek k=1
(j=m,i=l)

donde AEk
lm = AKgradŵm (x) gradŵl (x) dx es la matriz de conducción local y la expresión
Ek
n
#
Ek
Alm corresponde al proceso conocido como ensamble.
(j=m,i=l)
k=1
69
L a2 a
n+1
fi = f (t) wi (x) dΩh = f (t) wi (x) dA dx = Af (t) wi (x) dx + · · · + Af (t) wi (x) dx

Ωh 0 A a1 an
n
#
= Af (t) wi (x)|E1 dx + · · · + Af (t) wj (x)|En dx = Af (t) wi (x)|Ek dx
k=1 E
E1   En k
n
# n
#
=  Af (t) ŵl (x) dx = flEk
(i=l)
k=1 Ek k=1
(i=l)
n
#
donde flEk = Af (t) ŵl (x) dx es el vector de fuente interna local y la expresión Ek
fl
(i=l)
Ek k=1
corresponde al proceso conocido como ensamble.
n
#
qi = q̂ (t) γwi (x) dΓ2h = q̂ (t) γwi |Γ2h (x) dA = (A q̂ (t) γ ŵl (x))|Ek ∩Γ2h (i=l)
Γ2h A k=1
n
#
= qE
l
k ∩Γ2h
(i=l)
k=1
donde qEl
k ∩Γ2h
= A q̂ (t) γ ŵl (x)|Ek ∩Γ2h es el vector de cantidad de calor prescrita (cargas) local y la
# n
expresión qE
l
k ∩Γ2h
(i=l)
k=1
Sea Ek el elemento finito geométrico
el espacio de elemento finito local PEk es PEk = {p : Ek → R, p (x) = β 0 + β 1 x} , y el conjunto de

grados de libertad ΣEk es ΣEk = {φi ∈ L (PE , R) : φ1 (p) ∈ PEk , φ2 (p) ∈ PEk } .
De los resultados ya obtenidos previamente tenemos que el espacio de elemento finito local, PEk , está
dado por PEk = [w̃1 , w̃2 ] ,
b2 − x x − b1
w̃1 (x) = , x ∈ Ek , w̃2 (x) = , x ∈ Ek .
b2 − b1 b2 − b1
Consideremos l = 1, · · · , 2 y m = 1, · · · , 2, calculemos las matrices y vectores locales;

Ek
Mlm = cv ρAŵm (x) ŵl (x) dx, l = 1, 2, m = 1, 2.
Ek
70

b2 b2

Ek
M11 = cv ρAŵ1 (x) ŵ1 (x) dx MEk
= cv ρAŵ2 (x) ŵ2 (x) dx
22

b1 b1
b2 $ %$ % b2 $ %$ %
b2 − x b2 − x x − b1 x − b1
= cv ρA dx = cv ρA dx
b2 − b1 b2 − b1 b2 − b1 b2 − b1
b1 b1
$ % b2 $ % b2

cv ρA 2 cv ρA
= (b 2 − x) dx = (x − b1 )2 dx
(b2 − b1 ) 2 (b2 − b1 )2

$ % b 1 $ % b1
cv ρA 3 cv ρA
= 1
(b 2 − b1 )
= 1
(b2 − b1 )3
(b2 − b1 )2 3 (b2 − b1 ) 2 3
1
= 3 cv ρA (b2 − b1 ) = 13 cv ρA (b2 − b1 )

2
b b2
Ek
M12 = cv ρAŵ2 (x) ŵ1 (x) dx ME 21 =
k
cv ρAŵ1 (x) ŵ2 (x) dx

b1 b1
b2 $ %$ %
b2 $ %$ %
x − b1 b2 − x b2 − x x − b1
= cv ρA dx
b2 − b1 b2 − b1 = cv ρA b − b b2 − b1
dx
2 1
b1 b1
$ % b2 $ % b2
cv ρA cv ρA

= (x − b1 ) (b2 − x) dx = (b2 − x) (x − b1 ) dx
(b2 − b1 )2
(b 2 − b1 ) 2
$ % b1 $ % b1
cv ρA 1 3 1 2 1 3 1 2 cv ρA 1 3 1 2 1 3 1 2
= 2 b − 2 b1 b2 − 6 b1 + 2 b2 b1 =
6 2 2 b − 2 b1 b2 − 6 b1 + 2 b2 b1
6 2
$ (b 2 − b 1 ) % $ (b 2 − b1 ) %
cv ρA cv ρA
1 3 2 2 3
(b2 − 3b2 b1 + 3b2 b1 − b1 ) 1 3 2 2 3
=
(b − b ) 2 6 = (b − b )2 6 (b2 − 3b2 b1 + 3b2 b1 − b1 )
$ 2 1 % $ 2 1 %
cv ρA 3 cv ρA (b2 − b1 ) cv ρA
= 2
1
6
(b2 − b1 ) = = 2
1
6
(b2 − b1 )3 = 16 cv ρA (b2 − b1 )
(b2 − b1 ) 6 (b2 − b1 )
Finalmente la matriz de acumulación local está dada como:

Ek cv ρA (b2 − b1 ) 2.0 1.0
M =
6 1.0 2.0
long
si la longitud del elemento Ek es constante, b2 − b1 = entonces tenemos
n

Ek cv ρA (long) 2.0 1.0
M =
6n 1.0 2.0

AEk
lm = AKgradŵm (x) gradŵl (x) dx, l = 1, 2, m = 1, 2.
Ek
71

b2 b2
Ek Ek
A11 = AK grad ŵ1 (x) grad ŵ1 (x) dx A12 = AK grad ŵ2 (x) grad ŵ1 (x) dx

b1 b1
2
b $ % $ %
b2 $ % $ %
b2 − x b2 − x x − b1 b2 − x
= AK grad grad
dx = AK grad grad dx
b2 − b1 b2 − b1 b2 − b1 b2 − b1
b1 b1
b2 $ %$ %
b2 $ %$ %

−1 −1 1 −1
= Ak dx = Ak dx
b2 − b1 b2 − b1 b2 − b1 b2 − b1
b1 b1

$ % b2 $ % b2
Ak Ak
= 1 dx = −1 dx
(b2 − b1 )2 (b2 − b1 )2

$ % b1 $ % b1
Ak Ak Ak −Ak
= (b2 − b1 ) = = (− (b2 − b1 )) =
2
(b2 − b1 ) 2
(b2 − b1 )
(b2 − b1 ) (b2 − b1 )

2
b 2
b
Ek Ek
A21 = AK grad ŵ1 (x) grad ŵ2 (x) dx
A22 = AK grad ŵ2 (x) grad ŵ2 (x) dx

b1 b1
b2 $ % $ % b2 $ % $ %
b2 − x x − b1 x − b1 x − b1
= AK grad grad
dx = AK grad grad dx
b2 − b1 b2 − b1 b2 − b1 b2 − b1
b1 b1
b2 $ %$ % b2 $ %$ %

−1 1 1 1
= Ak dx = Ak dx
b2 − b1 b2 − b1 b2 − b1 b2 − b1
b1 b1

$ % b2 $ % b2
Ak Ak
= −1 dx = 1 dx
(b2 − b1 )2 (b2 − b1 )2

$ % b1 $ % b1
Ak −Ak Ak Ak
= (− (b2 − b1 )) = = (b2 − b1 ) =
2
(b2 − b1 ) 2
(b2 − b1 )
(b2 − b1 ) (b2 − b1 )
Finalmente la matriz de conductividad térmica local está dada como:

Ek Ak 1.0 −1.0
A =
(b2 − b1 ) −1.0 1.0
long
n

kAn 1.0 −1.0
AEk =
(long) −1.0 1.0

flEk = Af (t) ŵl (x) dx, l = 1, 2
Ek
72

b2 b2

Ek
f1 = Af (t) ŵ1 (x) dx f k = Af (t) ŵ (x) dx
E
2 2

b2 b1 $ % b1 $ %
b2
b2 − x x − b1
= Af (t) dx = Af (t) dx
b2 − b 1 b 2 − b 1
$ b1 % b2 $ b 1 % b2
Af (t) Af (t)
= (b2 − x) dx = (x − b1 ) dx

$ b2 − b1 % b1 $ b2 − b1 % b1
Af (t) 1 2
= Af (t)
1 2
= b − b2 b1 + 12 b21
2 2 2
b 2 − b 2 b 1 + 1 2
2
b 1
$ b2 − b1 % $ b2 − b1 %
Af (t) 1 2
= (b − 2b b + b 2
) = Af (t) 1 (b22 − 2b2 b1 + b21 )
2 2 2 1 1 2
$ b2 − b1 % $ b2 − b1 %
Af (t) 1
Af (t) (b2 − b1 ) Af (t) 1 Af (t) (b2 − b1 )
= 2
(b2 − b1 )2 = = 2
(b2 − b1 )2 =
b2 − b1 2 b2 − b1 2
Finalmente el vector fuente local está dado como:

Ek Af (t) (b2 − b1 ) 1.0
f = .
2 1.0
long
n

Ek Af (t) (long) 1.0
f = .
2n 1.0
qE k ∩Γ2h
= (A q̂ (t) γ ŵl (x))|Ek ∩Γ2h , l = 1, 2
l $ %
b2 − x A q̂ (t) si b1 ∈ Ek ∩ Γ2h
A q̂ (t) γ ŵ1 (x)|Ek ∩Γ2h = A q̂ (t) =
b2 − b1 Ek ∩Γ2h 0 si b2 ∈ Ek ∩ Γ2h
$ %
x − b1 0 si b1 ∈ Ek ∩ Γ2h
A q̂ (t) γ ŵ2 (x)|Ek ∩Γ2h = A q̂ (t) =
b2 − b1 Ek ∩Γ2h A q̂ (t) si b2 ∈ Ek ∩ Γ2h
Finalmente el vector de cantidad de calor (cargas) local está dado como:


 A q̂ (t)

 si b1 ∈ Ek ∩ Γ2h ,

 0

0
qEk ∩Γ2h = si Ek ∩ Γ2h = ∅. .

 0



 0
 si b2 ∈ Ek ∩ Γ2h .
A q̂ (t)
A partir de las matrices y vectores locales procedemos a realizar el ensamble de la matrices globales
long
para n = 10, con longitud de los elementos constante b2 − b1 = entonces tenemos .
n
73
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
M=  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

global j= , j=
local m
= 1, m = 2
cv ρA (long) 2.0 1.0 l=1 i=
MEk =
6n 1.0 2.0 l=2 i=
local global
74
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
A=  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

global j= , j=
local m
= 1, m = 2
kAn 1.0 −1.0 l=1 i=
AEk =
(long) −1.0 1.0 l=2 i=
local global
75
   

 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 

f= q=

 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


 
 
 


Ek Af (t) (long) 1.0 l=1 i=
f =
2n 1.0 l=2 i=
local global


 1.0 l=1 i=

 si b1 ∈ Ek ∩ Γ2h ,

 0 l=2 i=

0 l =1 i=
qEk ∩Γ2h = A q̂ (t) si Ek ∩ Γ2h = ∅.

 0 l=2 i=



 0 l=1 i=
 si b2 ∈ Ek ∩ Γ2h .
1.0 l=2 i=
local global
TAREA 3. Repetir el proceso de ensamble para las matrices y vectores del problema de conducción
de calor para la discretización dada:
76
Considere n = 10 y que las longitudes de los elementos es constante y la longitud total es long.

TAREA 4. Repetir el proceso de ensamble para las matrices y vectores del problema de conducción
de calor para la discretización dada:
ar+1 − ar
Considere n = 10 y que las longitudes para elementos adyacentes cumplen = 1.2 y la
ar+2 − ar+1
longitud total es long.
2.8. Conducción de Calor Bidimensional.
En el problema de conducción de calor evolutivo, consideremos el cuerpo Ω con las dimensiones

geométricas:
77
consideremos el elemento finito local:
Fig. 2.7. Elem. Fin. Geométrico.
con los valores de las propiedades físicas y dimensiones siguientes:


, ρ = 7850.
e ≡ espesor del cuerpo [m] , e = 0.01.
long = L ≡ altura y base del cuerpo [m] , long = 1.000.
n ≡ número de elementos finitos por lado del cuerpo en la discretización.
long
b = bx2 − bx1 ≡ base de los elementos finitos [m] , b = .
n.
long
h = by3 − by1 ≡ altura de los elementos finitos [m] , h = .
W n.
78
L L e
Mij = cv ρwj (x, y) wi (x, y) dΩh = cv ρwj (x, y) wi (x, y) dz dy dx
Ωh 0 0 0
L L
= ecv ρwj (x, y) wi (x, y) dy dx
0 0

= ecv ρwj (x, y) wi (x, y) dE1 + · · · + ecv ρwj (x, y) wi (x, y) dEn2
E1 En2
 
n2
# n2
#
= ecv ρ wj (x, y)|Ek wi (x, y)|Ek dEk =  ecv ρŵm (x, y) ŵl (x, y) dEk 
k=1 Ek k=1 Ek (j=m,i=l)
#n2
= MEk
lm
(j=m,i=l)
k=1

donde MEk
lm = ecv ρŵm (x, y) ŵl (x, y) dEk es la matriz de acumulación local y la expresión
Ek
n2
#
MEk
(j=m,i=l)
k=1
L L e
Aij = K grad wj (x, y) · grad wi (x, y) dΩh = K grad wj (x, y) · grad wi (x, y) dz dy dx
Ωh 0 0 0
L L
= eK grad wj (x, y) · grad wi (x, y) dy dx
0 0

= eK grad wj (x, y) · grad wi (x, y) dE1 + · · · + eK grad wj (x, y) · grad wi (x, y) dEn2
E1 En2
 
# n2
=  eK grad wj (x, y)| · grad wi (x, y)|Ek dEk 
Ek
k=1
Ek 
n2
# n
#
2

=  eKgradŵm (x, y) · gradŵl (x, y) dEk  = AEk
lm
(j=m,i=l)
k=1 Ek k=1
(j=m,i=l)

donde AEk
lm = eK grad ŵm (x, y)·grad ŵl (x, y) dEk es la matriz de conducción local y la expresión
Ek
n2
#
AEk
(j=m,i=l)
k=1
79
L L e L L
fi = f (t) wi (x, y) dΩh = f (t) wi (x, y) dz dy dx = ef (t) wi (x, y) dy dx
Ωh 0 0 0 0 0

= ef (t) wi (x, y) dE1 + · · · + ef (t) wi (x, y) dEn2
E1 En2
   
#
n2 n2
# n
#
2

=  ef (t) wi (x, y)| dEk  =  ef (t) ŵl (x, y) dEk  = flEk
Ek
(i=l)
k=1 Ek k=1 Ek k=1
(i=l)
n
# 2
donde flEk = ef (t) ŵl (x, y) dEk es el vector de fuente interna local y la expresión Ek
fl
(i=l)
Ek k=1
e
qi = q̂ (t) γwi (x, y) dΓ2h = q̂ (t) γwi (x, y) dz d (Ek ∩ Γ2h )
Γ Ek ∩Γ2h 0
2h
= eq̂ (t) γwi (x, y) d (E1 ∩ Γ2h ) + · · · + eq̂ (t) γwi (x, y) d (En2 ∩ Γ2h )
E1 ∩Γ2h En2 ∩Γ2h
 
n2
#
=  eq̂ (t) ( γwi (x, y))|Ek ∩Γ2h d (Ek ∩ Γ2h )
k=1
 Ek ∩Γ2h 
n2
# n
#
2

=  eq̂ (t) ( γ ŵl (x, y))|Ek ∩Γ2h d (Ek ∩ Γ2h ) = qE
l
k ∩Γ2h
(i=l)
k=1 Ek ∩Γ2h k=1
(i=l)

donde qE
l
k ∩Γ2h
= eq̂ (t) ( γ ŵl (x, y))|Ek ∩Γ2h d (Ek ∩ Γ2h ) es el vector de cantidad de calor pre-
Ek ∩Γ2h
n
# 2
scrita (cargas) local y la expresión qE

l
k ∩Γ2h
(i=l)
k=1
el espacio de elemento finito local PEk es PEk = {p : Ek → R, p (x, y) = β 0 + β 1 x + β 2 y + β 3 xy}

80
De los resultados ya obtenidos previamente tenemos que el espacio de elemento finito local, PEk , está
dado por PEk = [w̃1 , w̃2 , w̃3 , w̃4 ] ,
(bx2 by3 ) − (by3 ) x − (bx2 ) y + xy
w̃1 (x, y) = , (x, y) ∈ Ek ,
− (bx1 by3 ) + (by3 ) x + (bx1 ) y − xy
w̃2 (x, y) = , (x, y) ∈ Ek ,
− (bx2 by1 ) + (by1 ) x + (bx2 ) y − xy
w̃3 (x, y) = , (x, y) ∈ Ek ,
(bx1 by1 ) − (by1 ) x − (bx1 ) y + xy
w̃4 (x, y) = , (x, y) ∈ Ek ,
especializamos las funciones base de acuerdo a b1 = (0, 0) , b2 = (b, 0) , b3 = (0, h) y b4 = (b, h) ,
obteniéndose
bh − hx − by + xy hx − xy by − xy xy
w̃1 (x, y) = , w̃2 (x, y) = , w̃3 (x, y) = , w̃4 (x, y) = ,
bh bh bh bh
consideremos l = 1, · · · , 4 y m = 1, · · · , 4, calculemos las matrices y vectores locales;

Ek
Mlm = ecv ρŵm (x, y) ŵl (x, y) dEk , l = 1 · · · , 4, m = 1 · · · , 4.
Ek
h b h b $ %2
bh − hx − by + xy
MEk
11 = ecv ρŵ1 (x, y) ŵ1 (x, y) dx dy = ecv ρ dx dy
bh
0 0 0 0
h b
ecv ρ ecv ρ 1 3 3 ecv ρbh
= (bh − hx − by + xy)2 dx dy = bh =
(bh)2 (bh)2 9 9
0 0
h b h b $ %$ %
hx − xy bh − hx − by + xy
MEk
bh bh
0 0 0 0
h b
= (hx − xy) (bh − hx − by + xy) dx dy = bh =
(bh)2 (bh)2 18 18
0 0
h b h b $ %$ %
by − xy bh − hx − by + xy
MEk
bh bh
0 0 0 0
h b
= (by − xy) (bh − hx − by + xy) dx dy = bh =
(bh)2 (bh)2 18 18
0 0
h b h b $ %
xy bh − hx − by + xy
MEk
bh bh
0 0 0 0
h b
= (xy) (bh − hx − by + xy) dx dy = bh =
(bh)2 (bh)2 36 36
0 0
81
h b h b $ %$ %
hx − xy hx − xy
MEk
bh bh
0 0 0 0
h b
= (hx − xy) (hx − xy) dx dy = bh =
(bh)2 (bh)2 9 9
0 0
h b h b $ %$ %
by − xy hx − xy
MEk
bh bh
0 0 0 0
h b
= (by − xy) (hx − xy) dx dy = bh =
(bh)2 (bh)2 36 36
0 0
h b h b $ %
xy hx − xy
MEk
bh bh
0 0 0 0
h b
= (xy) (hx − xy) dx dy = bh =
(bh)2 (bh)2 18 18
0 0
h b h b $ %$ %
by − xy by − xy
MEk
bh bh
0 0 0 0
h b
= (by − xy) (by − xy) dx dy = bh =
(bh)2 (bh)2 9 9
0 0
h b h b $ %
xy by − xy
MEk
bh bh
0 0 0 0
h b
= (xy) (by − xy) dx dy = bh =
(bh)2 (bh)2 18 18
0 0
h b h b
xy xy
ME
44 =
k
ecv ρŵ4 (x, y) ŵ4 (x, y) dx dy = ecv ρ dx dy
bh bh
0 0 0 0
h b
= (xy) (xy) dx dy = bh =
(bh)2 (bh)2 9 9
0 0
 
4 2 2 1

e cv ρ b h  2 4 1 2 
MEk = 
36  2 1 4 2 
1 2 2 4
82
long
si b = h = simplificamos a:
n
 
4 2 2 1
2 
ecv ρ (long)  2 4 1 2 

MEk =  2
36n2 1 4 2 
1 2 2 4

AEk
lm = eK gradŵm (x, y) gradŵl (x, y) dEk , l = 1, · · · , 4, m = 1, · · · , 4.
Ek
h b
AE
11 =
k
eK gradŵ1 (x, y) · gradŵ1 (x, y) dx dy
0 0
h b $ % $ %
bh − hx − by + xy bh − hx − by + xy
= eK grad · grad dx dy
bh bh
0 0  T  
b
h
 −h + y 
    −h + y   h b $
h−yh−y b−xb−x
%
= ek bh bh
 −b + x   −b + x  dx dy = ek
bh bh
+
bh bh
dx dy
   
0 0 0 0
bh bh

h b
ek 2 2 ek 1 2 2
ek (h2 + b2 )
= (h − y) + (b − x) dx dy = (h + b ) bh =
(bh)2 (bh)2 3 3bh
0 0
h b
AE
12 =
k
0 0
h b % $ $ %
hx − xy bh − hx − by + xy
bh bh
0 0  T  
h b 
 h − y 
   −h + y 
 h b $ %
bh bh h − y −h + y −x −b + x
= ek −x   −b + x  dx dy = ek + dx dy

    bh bh bh bh
0 0 bh 0 0
bh
h b
ek 2 ek 1 2 2
ek (−2h2 + b2 )
= − (h − y) − x (−b + x) dx dy = (−2h + b ) bh =
(bh)2 (bh)2 6 6bh
0 0
h b
AE
13 =
k
0 0
h b $ % $ %
by − xy bh − hx − by + xy
bh bh
0 0
83
   
h b  −y T  −h + y  h b $ %
    −y −h + y b − x −b + x
= ek bh bh
 b − x   −b + x  dx dy = ek
bh bh
+
bh bh
dx dy
   
0 0 bh 0 0
bh

h b
ek 2 ek 1 2 2
ek (h2 − 2b2 )
= −y (−h + y) − (b − x) dx dy = (h − 2b ) bh =
(bh)2 (bh)2 6 6bh
0 0
h b
AE
14 =
k
0 0
h b xy $ %
bh bh
0 0
 y T  −h + y 
h b     
 h b $
y −h + y x −b + x
%
= ek bh bh
 x   −b + x  dx dy = ek
bh bh
+
bh bh
dx dy
 
0 0 bh bh 0 0

h b
ek ek 1 2 2
−ek (h2 + b2 )
= (y (−h + y) + x (−b + x)) dx dy = − (h + b ) bh =
(bh)2 (bh)2 6
6bh
0 0
h b
AE
22 =
k
0 0
h b %$ $ %
hx − xy hx − xy
bh bh
0 0  T  
h b 
 h − y 
   h − y 
 h b $ %
bh bh h − y h − y x x
= ek −x   −x  dx dy = ek + dx dy

    bh bh bh bh
0 0 bh bh 0 0
h b
ek 2 2
ek 1 2 2
ek (h2 + b2 )
= (h − y) + x dx dy = (h + b ) bh =
(bh)2 (bh)2 3 3bh
0 0
h b
AE
23 =
k
0 0
h b % $ $ %
by − xy hx − xy
= eK grad · gradŵ dx dy
bh bh
0 0    
b
h  −y T  h − y  h b $ %
    −y h − y b − x −x
= ek bh bh
 b − x   −x  dx dy = ek
bh bh
+
bh bh
dx dy
   
0 0 bh bh 0 0
84
h b
ek ek 1 2 2
−ek (h2 + b2 )
= (−y (h − y) − (b − x) x) dx dy = − (h + b ) bh =
(bh)2 (bh)2 6
6bh
0 0
h b
AE
24 =
k
0 0
h b b xy $ %
hx − xy
= · grad
eK grad dx dy
bh
0 bh
0 0
 y T  
b
h
    h−y   h b $
y h − y x −x
%
= ek bh bh dx dy = ek + dx dy
 x    −x   bh bh bh bh
0 0 bh bh 0 0
h b
ek 2 ek 1 2 2
ek (h2 − 2b2 )
= (y (h − y) − x ) dx dy = (h − 2b ) bh =
(bh)2 (bh)2 6 6bh
0 0
h b
AE
33 =
k
0 0
h b % $ $ %
by − xy by − xy
bh bh
0 0    
b
h  −y T  −y  h b $ %
    −y −y b − x b − x
= ek bh bh
 b − x   b − x  dx dy = ek
bh bh
+
bh bh
dx dy
   
0 0 bh bh 0 0

h b
ek 2 2 ek 1 2 2
ek (h2 + b2 )
= y + (b − x) dx dy = (h + b ) bh =
(bh)2 (bh)2 3 3bh
0 0
h b h b xy $ %
by − xy
AEk
34 = eK gradŵ4 (x, y) · gradŵ3 (x, y) dx dy = eK grad · grad dx dy
bh bh
 y T  −y
0 0  0 0
h b     
 h b $
y −y x b−x
%
= ek bh bh dx dy = ek + dx dy
 x   
b−x 
 bh bh bh bh
0 0 bh bh 0 0
h b
ek ek 1 ek (−2h2 + b2 )
= (−y 2 + x (b − x)) dx dy = (−2h2
+ b2
) bh =
(bh)2 (bh)2 6 6bh
0 0
h b h b xy
AE
44 =
k
eK gradŵ4 (x, y) · gradŵ4 (x, y) dx dy = eK grad · gradŵ4 (x, y) dx dy
bh
0 0 0 0
85
 y T  y 
h b    y y
bh bh x x
= ek dx dy = vek + dx dy
 x   x  bh bh bh bh
0 0 bh bh
h b
ek 2 2 ek 1 2 2
ek (h2 + b2 )
= (y + x ) dx dy = (h + b ) bh =
(bh)2 (bh)2 3 3bh
0 0
 
2 (h2 + b2 ) (−2h2 + b2 ) (h2 − 2b2 ) − (h2 + b2 )
ek 
 (−2h
2
+ b2 ) 2 (h2 + b2 ) − (h2 + b2 ) (h2 − 2b2 ) 
AEk =
6bh  (h2 − 2b2 ) − (h2 + b2 ) 2 (h2 + b2 ) (−2h2 + b2 ) 
− (h2 + b2 ) (h2 − 2b2 ) (−2h2 + b2 ) 2 (h2 + b2 )
long
n
 
4 −1 −1 −2
ek  
 −1 4 −2 −1 
AEk =
6  −1 −2 4 −1 
−2 −1 −1 4

flEk = ef (t) ŵl (x, y) dEk , l = 1, · · · , 4
Ek
h b h b $ %
f1Ek = ef (t) ŵ1 (x, y) dx dy = ef (t) dx dy
bh
0 0 0 0
h b
ef (t) ef (t) 1 2 2 ef (t) bh
= (bh − hx − by + xy) dx dy = hb =
bh bh 4 4
0 0
h b h b $ %
hx − xy
bh
0 0 0 0
h b
ef (t) ef (t) 1 2 2 ef (t) bh
= (hx − xy) dx dy = hb =
bh bh 4 4
0 0
h b h b $ %
by − xy
bh
0 0 0 0
h b
ef (t) ef (t) 1 2 2 ef (t) bh
= (by − xy) dx dy = hb =
bh bh 4 4
0 0
h b h b xy
bh
0 0 0 0
86
h b
ef (t) ef (t) 1 2 2 ef (t) bh
= (xy) dx dy = hb =
bh bh 4 4
0 0
 

 1.0 

ef (t) bh  1.0 
f Ek =
4 
 1.0 

 
1.0
long
n
 

 1.0 

ef (t) (long)2  1.0 
f Ek = .
4n2 
 1.0 

 
1.0

qE
l
k ∩Γ2h
= eq̂ (t) ( γ ŵl (x, y))|Ek ∩Γ2h d (Ek ∩ Γ2h ) , l = 1, · · · , 4
 Ekb∩Γ2h




 eq̂ (t) ( γ ŵ1 (x, y))|y=0 dx si b1 , b2 ∈ Ek ∩ Γ2h





 0

 h



 eq̂ (t) ( γ ŵ1 (x, y))|x=b dy si b2 , b4 ∈ Ek ∩ Γ2h


qE
1
k ∩Γ2h
= 0

 b



 eq̂ (t) ( γ ŵ1 (x, y))|y=h dx si b3 , b4 ∈ Ek ∩ Γ2h



 0



 h




 eq̂ (t) ( γ ŵ1 (x, y))|x=0 dy si b1 , b3 ∈ Ek ∩ Γ2h
0
 b 
 $ %  b

 bh − hx − by + xy 
 eq̂ (t)

 eq̂ (t) dx 
 (bh − hx) dx

 bh 
 bh

 y=0 


 0 
 0


 h $ % 
 h eq̂ (t) 1 2


 bh − hx − by + xy 

 eq̂ (t) 

 bh
 eq̂ (t) dy  0 dy  2

 bh 
 bh  bh

0
x=b
0 0
= $ % = =

 b 
 b 
 0

 bh − hx − by + xy 
 eq̂ (t) 




eq̂ (t)
bh dx 

 bh
0 dx  eq̂ (t) 1 bh2


 y=h 
 bh 2

 0 
 0

 h $ % 
 h

 bh − hx − by + xy 
 eq̂ (t)

 eq̂ (t) dy 
 (bh − by) dy

 bh 
 bh
x=0
0 0
87


 eq̂ (t)

 b si b1 , b2 ∈ Ek ∩ Γ2h
 2

0 si b2 , b4 ∈ Ek ∩ Γ2h
=

 0 si b3 , b4 ∈ Ek ∩ Γ2h


 eq̂ (t) h

si b1 , b3 ∈ Ek ∩ Γ2h
2   b
 b  $ %

 
 hx − xy


 eq̂ (t) ( γ ŵ2 (x, y))|y=0 dx si b1 , b2 ∈ Ek ∩ Γ2h 

 eq̂ (t) dx

 
 bh y=0

 0 
 0

 h 
 h $ %

 
 hx − xy

 


 eq̂ (t) ( γ ŵ2 (x, y))|x=b dy si b2 , b4 ∈ Ek ∩ Γ2h 
 eq̂ (t) dy
  bh x=b
Ek ∩Γ2h 0 0
q2 = = $ %

 b 
 b

 
 hx − xy

 eq̂ (t) ( γ ŵ2 (x, y))|y=h dx si b3 , b4 ∈ Ek ∩ Γ2h 
 eq̂ (t) dx

 
 bh

 0 
 0
y=h

 
 $ %

 h 
 h

 
 hx − xy


 eq̂ (t) ( γ ŵ2 (x, y))|x=0 dy si b1 , b3 ∈ Ek ∩ Γ2h 

 eq̂ (t) dy
bh x=0
0 0

 b

 eq̂ (t)

 x dx

 b



 0
 

 h
eq̂ (t) 1 2 eq̂ (t)

 eq̂ (t) 
 


 (h − y) dy 
 2
b 
 b si b1 , b2 ∈ Ek ∩ Γ2h

 h  eq̂b(t)



 2
eq̂ (t)
1 2
= 0
= h = h si b2 , b4 ∈ Ek ∩ Γ2h

 b 
 h 2

 2

 eq̂ (t) 
 0 
 0 si b3 , b4 ∈ Ek ∩ Γ2h

 0 dx 
 


 bh 0 0 si b1 , b3 ∈ Ek ∩ Γ2h

 0



 h

 eq̂ (t)

 0 dy
 bh
0 b  b
  $ %

 
 by − xy


 eq̂ (t) ( γ ŵ3 (x, y))|y=0 dx si b1 , b2 ∈ Ek ∩ Γ2h 

 eq̂ (t) dx

 
 bh y=0

 0 
 0

 h 
 h $ %

 
 by − xy

 


 eq̂ (t) ( γ ŵ3 (x, y))|x=b dy si b2 , b4 ∈ Ek ∩ Γ2h 
 eq̂ (t) dy
  bh x=b
qE
3
k ∩Γ2h
= 0
= 0
$ %

 b

 b

 
 by − xy

 eq̂ (t) ( γ ŵ3 (x, y))|y=h dx si b3 , b4 ∈ Ek ∩ Γ2h 
 eq̂ (t) dx

 
 bh

 0 
 0
y=h

 


 h 
 h $ %

 
 by − xy


 eq̂ (t) ( γ ŵ3 (x, y))|x=0 dy si b1 , b3 ∈ Ek ∩ Γ2h 

 eq̂ (t) dy
bh x=0
0 0
88

 b

 eq̂ (t)

 0 dx

 bh



 0
 

 h

 eq̂ (t) 
 0 
 0 si b1 , b2 ∈ Ek ∩ Γ2h

 0 dy 
 


 bh  0
  0
 si b2 , b4 ∈ Ek ∩ Γ2h
0 eq̂ (t) 1 2 eq̂ (t)
= = b = b si b3 , b4 ∈ Ek ∩ Γ2h

 b 
 b 2 
 2

 eq̂ (t) 
 



 b
(b − x) dx  eq̂ (t) 1 h2
  eq̂ (t) h si b1 , b3 ∈ Ek ∩ Γ2h


 h 2
2

 0

 h

 eq̂ (t)

 y dy

 h
0 b  b
  xy

 



 eq̂ (t) ( γ ŵ 4 (x, y))| y=0 dx si b 1 , b2 ∈ E k ∩ Γ 2h 

 eq̂ (t) dx

 
 bh y=0

 0 
 0

 h 
 h

 
 xy

 


 eq̂ (t) ( γ ŵ 4 (x, y))| x=b dy si b 2 , b4 ∈ E k ∩ Γ 2h 
 eq̂ (t) dy
  bh x=b
qE
4
k ∩Γ2h
= 0
= 0

 b

 b
xy

 


 eq̂ (t) ( γ ŵ 4 (x, y))| y=h dx si b 3 , b4 ∈ E k ∩ Γ 2h 
 eq̂ (t) dx

 
 bh y=h

 0 
 0

 


 h 
 h
xy

 



 eq̂ (t) ( γ ŵ 4 (x, y))| x=0 dy si b 1 , b3 ∈ E k ∩ Γ 2h 

 eq̂ (t) dy
bh x=0
0 0

 b

 eq̂ (t)

 0 dx

 bh



 0
 

 h

 eq̂ (t) 
 0 
 0 si b1 , b2 ∈ Ek ∩ Γ2h

 y dy 
 eq̂ (t) 1 2 
 eq̂ (t)

 h 
 

2
h h si b2 , b4 ∈ Ek ∩ Γ2h
= 0
= h = 2

 b 
 eq̂ (t) 1 2 
 eq̂ (t)

 eq̂ (t) 
 2
b 
 b si b3 , b4 ∈ Ek ∩ Γ2h

 x dx 
 b 
 2

 b 0 0 si b1 , b3 ∈ Ek ∩ Γ2h

 0



 h

 eq̂ (t)

 0 dy
 bh
0
   

 b 
 
 0 

Ek ∩Γ2h eq̂ (t)  b

Ek ∩Γ2h eq̂ (t) 
h

q = si b1 , b2 ∈ Ek ∩ Γ2h , q = si b2 , b4 ∈ Ek ∩ Γ2h ,
2   0 
 2   0 
    
0 h
89
   
 0
   
 h 

eq̂ (t)  0  eq̂ (t) 
0

qEk ∩Γ2h = si b3 , b4 ∈ Ek ∩ Γ2h , q Ek ∩Γ2h
= si b1 , b3 ∈ Ek ∩ Γ2h ,
2   b  2  h 
   
 

b 0
long
n
   
 1 
   0 
 

eq̂ (t) long  1  eq̂ (t) long 1

qEk ∩Γ2h = si b1 , b2 ∈ Ek ∩ Γ2h , qEk ∩Γ2h = si b2 , b4 ∈ Ek ∩ Γ2h ,
2n 
 0  2n  0 
   
  
0 1
   

 0 
 
 1 
eq̂ (t) long   eq̂ (t) long   
0 0
qEk ∩Γ2h = si b3 , b4 ∈ Ek ∩ Γ2h , qEk ∩Γ2h = si b1 , b3 ∈ Ek ∩ Γ2h .
2n 
 1  2n  1 
   
  
1 0
TAREA 5. Concrete el proceso de ensamble para las matrices y vectores del problema de conducción
de calor bidimensional para la discretización dada:
Considere que las fronteras lateral derecha y lateral izquierda tienen la condición de temperatura
prescrita y las fronteras superior está perfectamente aislada y la frontera inferior tiene un flujo
prescrito positivo.
90
2.9. Conducción de Calor Tridimensional.
En el problema de conducción de calor evolutivo, consideremos el cuerpo Ω con las dimensiones

geométricas:
consideremos el elemento finito geométrico local con las coordenadas locales:

b1 = (0, 0, 0) , b2 = (b, 0, 0) , b3 = (0, h, 0) , b4 = (0, 0, a) , b5 = (b, h, 0) , b6 = (0, h, a) , b7 = (b, 0, a) y
b8 = (b, h, a) ,
91
Fig. 2.10. Elemento Finito Geométrico.
Los valores de las propiedades físicas y dimensiones siguientes:


, ρ = 7850.
long = L ≡ altura, ancho y base del cuerpo [m] , long = 1.000.
n ≡ número de elementos finitos por arista del cuerpo en la discretización.
long
b ≡ base de los elementos finitos [m] , b = .
n.
long
h ≡ altura de los elementos finitos [m] , h = .
n.
long
a ≡ ancho de los elementos finitos [m] , a = .
W n.
L L L
Mij = cv ρwj (x, y, z) wi (x, y, z) dΩh = cv ρwj (x, y, z) wi (x, y, z) dz dy dx
Ωh 0 0 0

= cv ρwj (x, y, z) wi (x, y, z) dE1 + · · · + cv ρwj (x, y, z) wi (x, y, z) dEn3
E1 En3
#n3
= cv ρ wj (x, y, z)|Ek wi (x, y, z)|Ek dEk
k=1 Ek
 
#n3 n
#
3

=  cv ρŵm (x, y, z) ŵl (x, y, z) dEk  = Ek
Mlm
(j=m,i=l)
k=1 Ek k=1
(j=m,i=l)

donde MEk
lm = cv ρŵm (x, y, z) ŵl (x, y, z) dEk es la matriz de acumulación local y la expresión
Ek
92
n
#
3

MEk
(j=m,i=l)
k=1

Aij = K grad wj (x, y, z) grad wi (x, y, z) dΩh
Ωh
L L L
= K grad wj (x, y, z) grad wi (x, y, z) dz dy dx
0 0 0

= K grad wj (x, y, z) grad wi (x, y, z) dE1 + · · · + K grad wj (x, y, z) grad wi (x, y, z) dEn3
E1 En3
 
# n3
=  Kgrad wj (x, y, z)| grad wi (x, y, z)|Ek dEk 
Ek
k=1
Ek 
n3
# n
#
3

=  Kgradŵm (x, y, z) gradŵl (x, y, z) dEk  = Ek
Alm
(j=m,i=l)
k=1 Ek k=1
(j=m,i=l)

donde AEk
lm = Kgradŵm (x, y, z) gradŵl (x, y, z) dEk es la matriz de conducción local y la expre-
Ek
n
#
3
Ek
sión Alm corresponde al proceso conocido como ensamble.
(j=m,i=l)
k=1
L L L
fi = f (t) wi (x, y, z) dΩh = f (t) wi (x, y, z) dz dy dx
Ωh 0 0 0

= f (t) wi (x, y, z) dE1 + · · · + f (t) wi (x, y, z) dEn3
E1 En3
   
n3
# n3
# n
#
3

=  f (t) wi (x, y, z)|E dEk  =  f (t) ŵl (x, y, z) dEk  = flEk
k
(i=l)
k=1 Ek k=1 Ek k=1
(i=l)
n
# 3
donde flEk = f (t) ŵl (x, y, z) dEk es el vector de fuente interna local y la expresión flEk
(i=l)
Ek k=1

qi = q̂ (t) γwi (x, y, z) dΓ2h = q̂ (t) γwi (x, y, z) d (Ek ∩ Γ2h )
Γ Ek ∩Γ2h
2h
= q̂ (t) γwi (x, y, z) d (E1 ∩ Γ2h ) + · · · + q̂ (t) γwi (x, y, z) d (En3 ∩ Γ2h )
E1 ∩Γ2h En3 ∩Γ2h
93
 
n3
#
=  q̂ (t) ( γwi (x, y, z))|Ek ∩Γ2h d (Ek ∩ Γ2h )
k=1
Ek ∩Γ2h 
n3
# n
#
3

=  q̂ (t) ( γ ŵl (x, y, z))|Ek ∩Γ2h d (Ek ∩ Γ2h ) = Ek ∩Γ2h
ql
(i=l)
k=1 Ek ∩Γ2h k=1
(i=l)

donde qE
l
k ∩Γ2h
= q̂ (t) ( γ ŵl (x, y, z))|Ek ∩Γ2h d (Ek ∩ Γ2h ) es el vector de cantidad de calor pre-
Ek ∩Γ2h
n
# 3
Ek ∩Γ2h
scrita (cargas) local y la expresión ql corresponde al proceso conocido como ensamble.
(i=l)
k=1
el espacio de elemento finito local PEk es

PEk = {p : Ek → R, p (x, y, x) = β 0 + β 1 x + β 2 y + β 3 z + β 4 xy + β 5 yz + β 6 xz + β 7 xyz}
y el conjunto de grados de libertad ΣEk es
ΣEk = {φi ∈ L (PEk , R) : φ1 (p) ∈ PEk , φ2 (p) ∈ PEk , · · · , φ8 (p) ∈ PEk } .
Procedemos a obtener el espacio de elemento finito local, PEk = [w̃1 , w̃2 , w̃3 , w̃4 , w̃5 , w̃6 , w̃7 , w̃8 ] ,
PEk = {v : Ek → R | v (x, y, z) = α1 w̃1 (x, y, z) + α2 w̃2 (x, y, z) + · · · + α8 w̃8 (x, y, z)} ,
φi (p (xj , yj , zj )) = δ ij entonces w̃i (x, y, z) = φi (p) (x, y, z) , i = 1, · · · , 8.
para este elemento finito geométrico procederemos a calcular los coeficientes de las ocho funciones
base
 en un solo proceso:   

 φ 1−8 (p (0, 0, 0)) = 1, 0, 0, 0, 0, 0, 0, 0 1 0 0 0 0 0 0 0 
 β0 


   

 φ
 1−8 (p (b, 0, 0)) = 0, 1, 0, 0, 0, 0, 0, 0  1 b 0 0 0 0 0 0  β1 


  1 0  


 φ (p (0, h, 0)) = 0, 0, 1, 0, 0, 0, 0, 0  h 0 0 0 0 0  β2 

 1−8  1 0  
φ1−8 (p (0, 0, a)) = 0, 0, 0, 1, 0, 0, 0, 0 0 a 0 0 0 0 β3
⇒
 1 b

 =

 φ1−8 (p (b, h, 0)) = 0, 0, 0, 0, 1, 0, 0, 0  h 0 bh 0 0 0  β4 


  1 0  


 φ 1−8 (p (0, h, a)) = 0, 0, 0, 0, 0, 1, 0, 0  h a 0 ha 0 0  β5 


  1 b  


 φ1−8 (p (b, 0, a)) = 0, 0, 0, 0, 0, 0, 1, 0 0 a 0 0 ba 0 
 β6 

  
φ1−8 (p (b, h, a)) = 0, 0, 0, 0, 0, 0, 0, 1 1 b h a bh ha ba bha β7
94
               

 1 
 
 0 
 
 0  
 0 
 
 0 
 
 0 
 
 0   
 0 


 
 
 
 
 
 
 
 
 
 
 
 
 
 
 


 0 
 
 1 
 
 0  
 0 
 
 0 
 
 0 
 
 0 
 
 0 


 
 
 
 
 
 
 
 
 
 
 
 
 
 
 


 0   0   1   0   0   0   0   0 
              
              

0 0 0 1 0 0 0 0
, , , , , , ,

 0 
 
 0 
 
 0  
 0 
 
 1 
 
 0 
 
 0   
 0 


 
 
 
 
 
 
 
 
 
 
 
 
 
 
 


 0 
 
 0 
 
 0 
 
 0 
 
 0 
 
 1 
 
 0 
 
 0 


 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 
 0     0    0   0
   0   0   1   0 
                       

0 0 0 0 0 0 0 1
 −1  
1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
 1 b 0 0 0 0 0 0    −1 1
0 0 0 0 0 0 
  −1 b b 
 1 0 h 0 0 0 0 0    0 1
0 0 0 0 0 
  −1h h 
 1 0 0 a 0 0 0 0   0 0 1
0 0 0 0 
  =  a1 −1 −1 a 
 1 b h 0 bh 0 0 0   0 1
0 0 0 
   bh bh bh bh 
 1 0 h a 0 ha 0 0    1 −1
0 ha ha −1
0 1
0 0 
  ha ha 
 1 b 0 a 0 0 ba 0   1 −1 0 −1 0 0 1
0 
ba ba ba ba
−1 1 1 1 −1 −1 −1 1
1 b h a bh ha ba bha bha bha bha bha bha bha bha bha
para (x, y, z) ∈ Ek las funciones base se explicitan como:
1 1
w̃1 (x, y, z) = bha
(b − x) (h − y) (a − z) , w̃2 (x, y, z) = bha (h − y) (a − z) x,
1 1 1
w̃3 (x, y, z) = bha
(b − x) (a − z) y, w̃4 (x, y, z) = bha (b − x) (h − y) z, w̃5 (x, y, z) = bha
(a − z) xy,
1 1 1
w̃6 (x, y, z) = bha
(b − x) yz, w̃7 (x, y, z) = bha (h − y) xz, w̃8 (x, y, z) = bha xyz,
consideremos
l = 1, · · · , 8 y m = 1, · · · , 8, calculemos las matrices y vectores locales;
Ek
Mlm = cv ρŵm (x, y, z) ŵl (x, y, z) dEk , l = 1 · · · , 8, m = 1 · · · , 8.
Ek
a h b
ME
11 =
k
cv ρŵ1 (x, y, z) ŵ1 (x, y, z) dx dy dz
0 0 0
a h b
1
2
= cv ρ bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
cv ρ cv ρ 1 3 3 3 cv ρbha
= ((b − x) (h − y) (a − z))2 dx dy dz = bha =
(bha)2 (bha)2 27 27
0 0 0
a h b
ME
12 =
k
0 0 0
a h b
1
1

= cv ρ bha
(h − y) (a − z) x bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
95
a h b
= x (b − x) (h − y)2 (a − z)2 dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
13 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) (a − z) y bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
= (b − x)2 (h − y) y (a − z)2 dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
14 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) (h − y) z bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
= (b − x)2 (h − y)2 (a − z) z dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
15 =
k
0 0 0
a h b
1
1

= cv ρ bha
(a − z) xy bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
= (b − x) x (h − y) y (a − z)2 dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b
ME
16 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) yz bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
= (b − x)2 (h − y) y (a − z) z dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
96
a h b
MEk
17 = cv ρŵ7 (x, y, z) ŵ1 (x, y, z) dx dy dz
0 0 0
a h b
1
1

= cv ρ bha
(h − y) xz bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
= (b − x) x (h − y)2 (a − z) z dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b
ME
18 =
k
0 0 0
a h b
1
1

= cv ρ bha
xyz bha
(b − x) (h − y) (a − z) dx dy dz
0 0 0
a h b
= (b − x) x (h − y) y (a − z) z dx dy dz = hab =
(bha)2 (bha)2 216 216
0 0 0
a h b
ME
22 =
k
0 0 0
a h b
1
2
= cv ρ bha
(h − y) (a − z) x dx dy dz
0 0 0
a h b
= ((h − y) (a − z) x)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
a h b
ME
23 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) (a − z) y bha
0 0 0
a h b
= ((b − x) (a − z) y) ((h − y) (a − z) x) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b
ME
24 =
k
0 0 0
97
a h b
1
1

= cv ρ bha
(b − x) (h − y) z bha
0 0 0
a h b
= ((b − x) (h − y) z) ((h − y) (a − z) x) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b
ME
25 =
k
0 0 0
a h b
1
1

= cv ρ bha
(a − z) xy bha
0 0 0
a h b
= ((a − z) xy) ((h − y) (a − z) x) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
26 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) yz bha
0 0 0
a h b
= ((b − x) yz) ((h − y) (a − z) x) dx dy dz = hab =
(bha)2 (bha)2 216 216
0 0 0
a h b
ME
27 =
k
0 0 0
a h b
1
1

= cv ρ bha
(h − y) xz bha
0 0 0
a h b
= ((h − y) xz) ((h − y) (a − z) x) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
28 =
k
0 0 0
a h b
1
1

= cv ρ bha
xyz bha
0 0 0
98
a h b
= (xyz) ((h − y) (a − z) x) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b a h b
2
ME
33 =
k
cv ρŵ3 (x, y, z) ŵ3 (x, y, z) dx dy dz = cv ρ 1
bha
(b − x) (a − z) y dx dy dz
0 0 0 0 0 0
a h b
= ((b − x) (a − z) y)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
a h b
ME
34 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) (h − y) z bha
0 0 0
a h b
= ((b − x) (h − y) z) ((b − x) (a − z) y) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b
ME
35 =
k
0 0 0
a h b
1
1

= cv ρ bha
(a − z) xy bha
0 0 0
a h b
= ((a − z) xy) ((b − x) (a − z) y) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
36 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) yz bha
0 0 0
a h b
= ((b − x) yz) ((b − x) (a − z) y) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
37 =
k
0 0 0
99
a h b
1
1

= cv ρ bha
(h − y) xz bha
0 0 0
a h b b
= ((h − y) xz) ((b − x) (a − z) y) dx dy dz = hab =
(bha)2 0 (bha)2 216 216
0 0 0
a h b
MEk
0 0 0
a h b
1
1

= cv ρ bha
xyz bha
0 0 0
a h b
= (xyz) ((b − x) (a − z) y) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b a h b
2
MEk
44 = cv ρŵ4 (x, y, z) ŵ4 (x, y, z) dx dy dz = cv ρ 1
bha
(b − x) (h − y) z dx dy dz
0 0 0 0 0 0
a h b
= ((b − x) (h − y) z)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
a h b
ME
45 =
k
0 0 0
a h b
1
1

= cv ρ bha
(a − z) xy bha
0 0 0
a h b
= ((a − z) xy) ((b − x) (h − y) z) dx dy dz = hab =
(bha)2 (bha)2 216 216
0 0 0
a h b
ME
46 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) yz bha
0 0 0
a h b
= ((b − x) yz) ((b − x) (h − y) z) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
100
a h b
MEk
0 0 0
a h b
1

= cv ρ bha
(h − y) xz ((b − x) (h − y) z) dx dy dz
0 0 0
a h b
= ((h − y) xz) ((b − x) (h − y) z) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b
ME
48 =
k
0 0 0
a h b
1
1

= cv ρ bha
xyz bha
0 0 0
a h b
= (xyz) ((b − x) (h − y) z) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b a h b
2
MEk
bha
(a − z) xy dx dy dz
0 0 0 0 0 0
a h b
= ((a − z) xy)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
a h b
ME
56 =
k
0 0 0
a h b
1
1

= cv ρ bha
(b − x) yz bha
0 0 0
a h b
= ((b − x) yz) ((a − z) xy) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b
ME
57 =
k
0 0 0
a h b
1
1

= cv ρ bha
(h − y) xz bha
0 0 0
101
a h b
= ((h − y) xz) ((a − z) xy) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b a h b

ME
58 =
k
bha
xyz 1
bha
0 0 0 0 0 0
a h b
= (xyz) ((a − z) xy) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b a h b
2
MEk
bha
(b − x) yz dx dy dz
0 0 0 0 0 0
a h b
= ((b − x) yz)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
a h b
ME
67 =
k
0 0 0
a h b
1
1

= cv ρ bha
(h − y) xz bha
0 0 0
a h b
= ((h − y) xz) ((b − x) yz) dx dy dz = hab =
(bha)2 (bha)2 108 108
0 0 0
a h b a h b

ME
68 =
k
bha
xyz 1
bha
0 0 0 0 0 0
a h b
= (xyz) ((b − x) yz) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b a h b
2
MEk
bha
(h − y) xz dx dy dz
0 0 0 0 0 0
a h b
= ((h − y) xz)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
a h b a h b

MEk
bha
xyz 1
bha
(h − y) xz dx dy dz
0 0 0 0 0 0
102
a h b
= (xyz) ((h − y) xz) dx dy dz = hab =
(bha)2 (bha)2 54 54
0 0 0
a h b a h b
2
ME
88 =
k
bha
xyz dx dy dz
0 0 0 0 0 0
a h b
= (xyz)2 dx dy dz = hab =
(bha)2 (bha)2 27 27
0 0 0
 
8 4 4 4 2 2 2 1
 4 8 2 2 4 1 4 2 
 
 4 2 8 2 4 4 1 2 
 
Ek cv ρbha 
 4 2 2 8 1 4 4 2 

M =
216   2 4 4 1 8 2 2 4 

 2 1 4 4 2 8 2 4 
 
 2 4 1 4 2 2 8 4 
1 2 2 2 4 4 4 8
long
si b = h = a = n
tenemos:
 
8 4 4 4 2 2 2 1
 4 8 2 2 4 1 4 2 
 
 4 2 8 2 4 4 1 2 
3  
cv ρ (long)  4 2 2 8 1 4 4 2 

MEk =
216n3   2 4 4 1 8 2 2 4 

 2 1 4 4 2 8 2 4 
 
 2 4 1 4 2 2 8 4 
1 2 2 2 4 4 4 8

AEk
lm = Kgrad ŵm (x, y, z) · grad ŵl (x, y, z) dEk , l = 1, · · · , 8, m = 1, · · · , 8.
Ek
a h b
AEk
11 = K gradŵ1 (x, y, z) · gradŵ1 (x, y, z) dx dy dz
0 0 0
a h b $ % $ %
(b − x) (h − y) (a − z) (b − x) (h − y) (a − z)
= kI grad · grad dx dy dz
bha bha
0 0 0
 T  
a h b  − (h − y) (a − z)   − (h − y) (a − z) 
k
= − (b − x) (a − z) − (b − x) (a − z) dx dy dz
(bha)2 
− (b − x) (h − y)
 
− (b − x) (h − y)

0 0 0
103
a h b
k 2 2 2 2 2 2
= (h − y) (a − z) + (b − x) (a − z) + (b − x) (h − y) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 (bh a + b ha + b h a) = (h2 a2 + b2 a2 + b2 h2 )
(bha) 9 9bha
a h b
AE
12 =
k
K gradŵ2 (x, y, z) · gradŵ1 (x, y, z) dx dy dz
0 0 0
a h b $ % $ %
x (h − y) (a − z) (b − x) (h − y) (a − z)
bha bha
0 0 0
   
a h b  (h − y) (a − z) T  − (h − y) (a − z) 
k
= −x (a − z) − (b − x) (a − z) dx dy dz
(bha)2 
−x (h − y)
 
− (b − x) (h − y)

0 0 0
a h b
k
= 2 − (h − y)2 (a − z)2 + x (b − x) (a − z)2 + x (b − x) (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 18 (−2bh a + b ha + b h a) = (−2h2 a2 + b2 a2 + b2 h2 )
(bha) 18bha
a h b
AE13 =
k
0 0 0
a h b $ % $ %
(b − x) y (a − z) (b − x) (h − y) (a − z)
bha bha
0 0 0    
a h b  −y (a − z) T  − (h − y) (a − z) 
k
= (b − x) (a − z) − (b − x) (a − z) dx dy dz
(bha)2 
− (b − x) y
 
− (b − x) (h − y)

0 0 0
a h b
k 2 2 2 2
= y (h − y) (a − z) − (b − x) (a − z) + (b − x) y (h − y) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 18 (b h a + h a b − 2b ha ) = (b2 h2 + h2 a2 − 2b2 a2 )
(bha) 18bha
a h b
AE
14 =
k
0 0 0
a h b $ % $ %
(b − x) (h − y) z (b − x) (h − y) (a − z)
bha bha
0 0 0
104
 T  
a h b  − (h − y) z   − (h − y) (a − z) 
k
= − (b − x) z − (b − x) (a − z) dx dy dz
(bha)2 
(b − x) (h − y)
 
− (b − x) (h − y)

0 0 0
a h b
k
= 2 (h − y)2 z (a − z) + (b − x)2 z (a − z) − (b − x)2 (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 18 (−2b h a + b ha + h a b) = (−2b2 h2 + b2 a2 + h2 a2 )
(bha) 18bha
a h b
Ek
A15 = K gradŵ5 (x, y, z) · gradŵ1 (x, y, z) dx dy dz
0 0 0
a h b % $ $ %
xy (a − z) (b − x) (h − y) (a − z)
bha bha
0 0 0    
a h b  y (a − z) T  − (h − y) (a − z) 
k
= x (a − z) − (b − x) (a − z) dx dy dz
(bha)2 
−xy
 
− (b − x) (h − y)

0 0 0
a h b
k
= 2 −y (h − y) (a − z)2 − x (b − x) (a − z)2 + x (b − x) y (h − y) dx dy dz
(bha)
0 0 0
k 1 2 2 2 2 2 2 k
= 2 36 (b h − 2h a − 2b a ) = (b2 h2 − 2h2 a2 − 2b2 a2 )
(bha) 36bha
a h b
Ek
0 0 0
a h b % $ $ %
(b − x) yz (b − x) (h − y) (a − z)
bha bha
0 0 0    
a h b  −yz T  − (h − y) (a − z) 
k
= (b − x) z − (b − x) (a − z) dx dy dz
(bha)2 
(b − x) y
 
− (b − x) (h − y)

0 0 0
a h b
k
= 2 y (h − y) z (a − z) − (b − x)2 z (a − z) − (b − x)2 y (h − y) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−2b h a − 2b ha + h a b) = (−2b2 h2 − 2b2 a2 + h2 a2 )
(bha) 36bha
a h b
Ek
0 0 0
105
a h b $ % $ %
x (h − y) z (b − x) (h − y) (a − z)
bha bha
0 0 0
   
a h b  (h − y) z T  − (h − y) (a − z) 
k
= −xz − (b − x) (a − z) dx dy dz
(bha)2 
x (h − y)
 
− (b − x) (h − y)

0 0 0
a h b
k 2 2
= − (h − y) z (a − z) + x (b − x) z (a − z) − x (b − x) (h − y) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−2b h a + b ha − 2h a b) = (−2b2 h2 + b2 a2 . − 2h2 a2 )
(bha) 36bha
a h b
AE
18 =
k
0 0 0
a h b xyz $ %
(b − x) (h − y) (a − z)
bha bha
0 0 0    
a h b  yz T  − (h − y) (a − z) 
k
= xz − (b − x) (a − z) dx dy dz
(bha)2 
xy
 
− (b − x) (h − y)

0 0 0
a h b
k
= (−y (h − y) z (a − z) − x (b − x) z (a − z) − x (b − x) y (h − y)) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−b h a − b ha − h a b) = (−b2 h2 − b2 a2 . − h2 a2 )
(bha) 36bha
a h b
Ek
a h 0 b 0 0 $ % $ %
x (h − y) (a − z) x (h − y) (a − z)
0 0 0  bha   bha 
a h b  (h − y) (a − z) T  (h − y) (a − z) 
k
= −x (a − z) −x (a − z) dx dy dz
(bha)2 0 0 0  −x (h − y)   −x (h − y) 
a h b
k
= 2 (h − y)2 (a − z)2 + x2 (a − z)2 + x2 (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 9 (b h a + h a b + b ha ) = (b2 h2 + h2 a2 + b2 a2 )
(bha) 9bha
a h b
AE
23 =
k
0 0 0
106
a h b $ % $ %
(b − x) y (a − z) x (h − y) (a − z)
bha bha
0 0 0
   
a h b  −y (a − z) T  (h − y) (a − z) 
k
= (b − x) (a − z) −x (a − z) dx dy dz
(bha)2 
− (b − x) y
 
−x (h − y)

0 0 0
a h b
k 2 2
= −y (h − y) (a − z) − x (b − x) (a − z) + x (b − x) y (h − y) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (b h a − 2h a b − 2b ha ) = (b2 h2 − 2h2 a2 − 2b2 a2 )
(bha) 36bha
a h b
AE
24 =
k
0 0 0
a h b $ % $ %
(b − x) (h − y) z x (h − y) (a − z)
bha bha
0 0 0    
a h b  − (h − y) z T  (h − y) (a − z) 
k
= − (b − x) z −x (a − z) dx dy dz
(bha)2 
(b − x) (h − y)
 
−x (h − y)

0 0 0
a h b
k
= 2 − (h − y)2 z (a − z) + x (b − x) z (a − z) − x (b − x) (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−2h ba + hb a − 2h b a) = (−2h2 a2 + b2 a2 − 2h2 b2 )
(bha) 36bha
a h b
Ek
0 0 0
a h b $ % $ %
xy (a − z) x (h − y) (a − z)
bha bha
0 0 0    
a h b  y (a − z) T  (h − y) (a − z) 
k
= x (a − z) −x (a − z) dx dy dz
(bha)2 
−xy
 
−x (h − y)

0 0 0
a h b
k
= 2 y (h − y) (a − z)2 − x2 (a − z)2 + x2 y (h − y) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 18 (h ba − 2hb a + h b a) = (h2 a2 − 2b2 a2 + h2 b2 )
(bha) 18bha
107
a h b
AEk
0 0 0
a h b % $ $ %
(b − x) yz x (h − y) (a − z)
bha bha
0 0 0    
a h b  −yz T  (h − y) (a − z) 
k
= (b − x) z −x (a − z) dx dy dz
(bha)2 
(b − x) y
 
−x (h − y)

0 0 0
a h b
k
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−h ba − h b a − hb a ) = (−h2 a2 − b2 a2 − h2 b2 )
(bha) 36bha
a h b
Ek
0 0 0
a h b % $ $ %
x (h − y) z x (h − y) (a − z)
bha bha
0 0 0    
a h b  (h − y) z T  (h − y) (a − z) 
k
= −xz −x (a − z) dx dy dz
(bha)2 
x (h − y)
 
−x (h − y)

0 0 0
a h b
k
= 2 (h − y)2 z (a − z) + x2 z (a − z) − x2 (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 18 (h ba + hb a − 2h b a) = (h2 a2 + b2 a2 − 2h2 b2 )
(bha) 18bha
a h b
AE
28 =
k
0 0 0
a h b xyz $ %
x (h − y) (a − z)
bha bha
0 0 0
   
a h b  yz T  (h − y) (a − z) 
k
= xz −x (a − z) dx dy dz
(bha)2 
xy
 
−x (h − y)

0 0 0
a h b
k
= 2 (y (h − y) z (a − z) − x2 z (a − z) − x2 y (h − y)) dx dy dz
(bha)
0 0 0
108
k 1 3 3 3 3 3 3 k
= 2 36 (h ba − 2h b a − 2hb a ) = (h2 a2 − 2b2 a2 − 2h2 b2 )
(bha) 36bha
a h b
AE
33 =
k
0 0 0
a h b $ % $ %
(b − x) y (a − z) (b − x) y (a − z)
bha bha
0 0 0    
a h b  −y (a − z) T  −y (a − z) 
k
= (b − x) (a − z) (b − x) (a − z) dx dy dz
(bha)2 
− (b − x) y
 
− (b − x) y

0 0 0
a h b
k 2
= 2 y (a − z)2 + (b − x)2 (a − z)2 + (b − x)2 y 2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 9 (h ba + hb a + h b a) = (h2 a2 + b2 a2 + h2 b2 )
(bha) 9bha
a h b
Ek
0 0 0
a h b $ % $ %
(b − x) (h − y) z (b − x) y (a − z)
bha bha
0 0 0    
a h b  − (h − y) z T  −y (a − z) 
k
= − (b − x) z (b − x) (a − z) dx dy dz
(bha)2 
(b − x) (h − y)
 
− (b − x) y

0 0 0
a h b
k
= 2 y (h − y) z (a − z) − (b − x)2 z (a − z) − (b − x)2 y (h − y) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (h ba − 2h b a − 2hb a ) = (h2 a2 − 2b2 a2 − 2h2 b2 )
(bha) 36bha
a h b
AE
35 =
k
0 0 0
a h b $ % $ %
xy (a − z)
(b − x) y (a − z)
bha bha
0 0 0
   
a h b  y (a − z) T  −y (a − z) 
k
= x (a − z) (b − x) (a − z) dx dy dz
(bha)2 
−xy
 
− (b − x) y

0 0 0
109
a h b
k 2 2 2 2

= −y (a − z) + x (b − x) (a − z) + x (b − x) y dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 18 (−2h ba + hb a + h b a) = (−2h2 a2 + b2 a2 + h2 b2 )
(bha) 18bha
a h b
AE
36 =
k
0 0 0
a h b %
$ $ %
(b − x) yz (b − x) y (a − z)
bha bha
0 0 0
   
a h b  −yz T  −y (a − z) 
k
= (b − x) z (b − x) (a − z) dx dy dz
(bha)2 
(b − x) y
 
− (b − x) y

0 0 0
a h b
k 2 2 2 2
= y z (a − z) + (b − x) z (a − z) − (b − x) y dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 18bha
a h b
AE37 =
k
0 0 0
a h b %$ $ %
x (h − y) z (b − x) y (a − z)
bha bha
0 0 0    
a h b  (h − y) z T  −y (a − z) 
k
= −xz (b − x) (a − z) dx dy dz
(bha)2 
x (h − y)
 
− (b − x) y

0 0 0
a h b
k
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 36bha
a h b
Ek
0 0 0
a h b xyz $ %
(b − x) y (a − z)
bha bha
0 0 0
110
   
a h b  yz T  −y (a − z) 
k
= xz (b − x) (a − z) dx dy dz
(bha)2 
xy
 
− (b − x) y

0 0 0
a h b
k
= 2 (−y 2 z (a − z) + x (b − x) z (a − z) − x (b − x) y 2 ) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 36bha
a h b
Ek
0 0 0
a h b $ % $ %
(b − x) (h − y) z (b − x) (h − y) z
bha bha
0 0 0    
a h b  − (h − y) z T  − (h − y) z 
k
= − (b − x) z − (b − x) z dx dy dz
(bha)2 
(b − x) (h − y)
 
(b − x) (h − y)

0 0 0
a h b
k
= 2 (h − y)2 z 2 + (b − x)2 z 2 + (b − x)2 (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 9bha
a h b
Ek
0 0 0
a h b $ % $ %
xy (a − z) (b − x) (h − y) z
bha bha
0 0 0    
a h b  y (a − z) T  − (h − y) z 
k
= x (a − z) − (b − x) z dx dy dz
(bha)2 
−xy
 
(b − x) (h − y)

0 0 0
a h b
k
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 36bha
a h b
AE
46 =
k
0 0 0
111
a h b $ % $ %
(b − x) yz (b − x) (h − y) z
bha bha
0 0 0
   
a h b  −yz T  − (h − y) z 
k
= (b − x) z − (b − x) z dx dy dz
(bha)2 
(b − x) y
 
(b − x) (h − y)

0 0 0
a h b
k 2 2 2 2
= y (h − y) z − (b − x) z + (b − x) y (h − y) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 18bha
a h b
AE
47 =
k
0 0 0
a h b $ % $ %
x (h − y) z (b − x) (h − y) z
bha bha
0 0 0    
a h b  (h − y) z T  − (h − y) z 
k
= −xz − (b − x) z dx dy dz
(bha)2 
x (h − y)
 
(b − x) (h − y)

0 0 0
a h b
k
= 2 − (h − y)2 z 2 + x (b − x) z 2 + x (b − x) (h − y)2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 18bha
a h b
Ek
0 0 0
a h b xyz $ %
(b − x) (h − y) z
bha bha
0 0 0    
a h b  yz T  − (h − y) z 
k
= xz − (b − x) z dx dy dz
(bha)2 
xy
 
(b − x) (h − y)

0 0 0
a h b
k
= 2 (−y (h − y) z 2 − x (b − x) z 2 + x (b − x) y (h − y)) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−2h ba − 2hb a + h b a) = (−2h2 a2 − 2b2 a2 + h2 b2 )
(bha) 36bha
112
a h b
AEk
0 0 0
a h b % $ $ %
xy (a − z) xy (a − z)
bha bha
0 0 0    
a h b  y (a − z) T  y (a − z) 
k
= x (a − z) x (a − z) dx dy dz
(bha)2 
−xy
 
−xy

0 0 0
a h b
k 2
= 2 y (a − z)2 + x2 (a − z)2 + x2 y 2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 9bha
a h b
Ek
0 0 0
a h b % $ $ %
(b − x) yz xy (a − z)
bha bha
0 0 0    
a h b  −yz T  y (a − z) 
k
= (b − x) z x (a − z) dx dy dz
(bha)2 
(b − x) y
 
−xy

0 0 0
a h b
k
= (−y 2 z (a − z) + x (b − x) z (a − z) − x (b − x) y 2 ) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 36bha
a h b
AE
57 =
k
0 0 0
a h b % $ $ %
x (h − y) z xy (a − z)
bha bha
0 0 0
   
a h b  (h − y) z T  y (a − z) 
k
= −xz x (a − z) dx dy dz
(bha)2 
x (h − y)
 
−xy

0 0 0
a h b
k
= 2 (y (h − y) z (a − z) − x2 z (a − z) − x2 y (h − y)) dx dy dz
(bha)
0 0 0
113
k 1 3 3 3 3 3 3 k
= 2 36 (h ba − 2h b a − 2hb a ) = (h2 a2 − 2h2 b2 − 2b2 a2 )
(bha) 36bha
a h b
AE
58 =
k
0 0 0
a h b xyz $ %
xy (a − z)
bha bha
0 0 0    
a h b  yz T  y (a − z) 
k
= xz x (a − z) dx dy dz
(bha)2 
xy
 
−xy

0 0 0
a h b
k
= 2 (y 2 z (a − z) + x2 z (a − z) − x2 y 2 ) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 18bha
a h b
Ek
0 0 0
a h b $
% $ %
(b − x) yz (b − x) yz
bha bha
0 0 0    
a h b  −yz T  −yz 
k
= (b − x) z (b − x) z dx dy dz
(bha)2 
(b − x) y
 
(b − x) y

0 0 0
a h b
k 2 2
= 2 y z + (b − x)2 z 2 + (b − x)2 y 2 dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 9bha
a h b
Ek
0 0 0
a h b % $ $ %
x (h − y) z (b − x) yz
bha bha
0 0 0
   
a h b  (h − y) z T  −yz 
k
= −xz (b − x) z dx dy dz
(bha)2 
x (h − y)
 
(b − x) y

0 0 0
114
a h b
k
= (−y (h − y) z 2 − x (b − x) z 2 + x (b − x) y (h − y)) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
= 2 36 (−2h ba − 2hb a + h b a) = (−2h2 a2 − 2b2 a2 + h2 b2 )
(bha) 36bha
a h b
AE
68 =
k
0 0 0
a h b xyz $ %
(b − x) yz
bha bha
0 0 0
   
a h b  yz T  −yz 
k
= xz (b − x) z dx dy dz
(bha)2 
xy
 
(b − x) y

0 0 0
a h b
k
= (−y 2 z 2 + x (b − x) z 2 + x (b − x) y 2 ) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 18bha
a h b
AE77 =
k
0 0 0
a h b $
% $ %
x (h − y) z x (h − y) z
bha bha
0 0 0    
a h b  (h − y) z T  (h − y) z 
k
= −xz −xz dx dy dz
(bha)2 
x (h − y)
 
x (h − y)

0 0 0
a h b
k 2 2 2 2 2 2
= (h − y) z + x z + x (h − y) dx dy dz
(bha)2
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 9bha
a h b
Ek
0 0 0
a h b xyz $ %
x (h − y) z
bha bha
0 0 0
115
   
a h b  yz T  (h − y) z 
k
= xz −xz dx dy dz
(bha)2 
xy
 
x (h − y)

0 0 0
a h b
k
= 2 (y (h − y) z 2 − x2 z 2 + x2 y (h − y)) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 18bha
a h b
Ek
0 0 0
   
a h b xyz xyz a h b  yz T  yz 
k
= kI grad · grad dx dy dz = xz xz dx dy dz
bha bha (bha)2 
xy
 
xy

0 0 0 0 0 0
a h b
k
= 2 (y 2 z 2 + x2 z 2 + x2 y 2 ) dx dy dz
(bha)
0 0 0
k 1 3 3 3 3 3 3 k
(bha) 9bha

Ek k [Ais ] [Ads ]
A =
36bha [Aii ] [Adi ]
donde:       2 2    
h2 a2 −2h2 a2 bh −2b2 h2
 4  +b2 a2  2  +b2 a2  2  +h2 a2  2  +b2 a2  
 
 2 2 2 2 2 2
a   +h2 a2  
  +b 2h 2   +b2 h2   −2b 
 −2h a bh b2 h2 −2h2 a2 
 
 2  +b2 a2  4  +h2 a2   −2h2 a2   +b2 a2  
 
 2 2 2 2
−2b2 a2   −2h2 b2  
[Ais ] =   +b2 h2   +b a   
 bh b2 h2 h2 a2 h2 a2 
 
 2  +h2 a2   −2h2 a2  4  +b2 a2   −2b2 a2  
 
 2 2 2 2 2 2
−2h2 b2  
  −2b2 a2   −2b 2a 2   +h b   
 −2b h −2h a h2 a2 h2 a2 
 
 2  +b a   +b a   −2b a  4  +b2 a2 
2 2 2 2 2 2 
+h2 a2 −2h2 b2 −2h2 b2 +h2 b2
116
         
b2 h2 −2b2 h2 −2b2 h2 −b2 h2
  −2h2 a2   −2b2 a2   +b2 a2 .   −b2 a2 .  
 
 2 2 2 2 2 2 2 2 
  −2b2 a2   +h2a2   −2h2 a2   −h2 a2  
 ha −h a ha ha 
 
 2  −2b a   −b a  2  +b a   −2b a  
2 2 2 2 2 2 2 2
 
 +h 2 2
b −h 2 2
b −2h 2 2
b −2h 2 2
b 
[Ads ] = 

 2 2
  2 2
  2 2
  2 2
 

 −2h a ha −h a −2h a 
 2  +b2 a2  2  +b2 a2   −b2 a2   +b2 a2  
 
 +h 2 2
b −2h 2 2
b −h 2 2
b −2h 2 2
b 
         
 −h a 2 2
ha2 2
−2h a 2 2
−2h a 2 2 
 
  −b2 a2  2  −2b2 a2  2  +b2 a2   −2b2 a2  
−h2 b2  +h2 b2 +h2 b2  +h2 b2
    
b2 h 2 h2 a2 −2h2 a2 −h2 a2
  −2h2 a2  2  −2b2 a2  2  +b2 a2   −b2 a2  
 
 −2b 2 2
a +h 2 2
b +h 2 2
b −h 2 2
b 
         
 2
−2b h 2 2
−h a 2 2
ha 2 2
ha 2 
 
  −2b2 a2   −b2 a2  2  +b2 a2  2  −2b2 a2  
 
 +h 2 2
a −h 2 2
b −2h 2 2
b +h 2 2
b 
[Aii ] = 

 2 2
  2 2
  2 2
  2 2
 

 −2b h ha −h a −2h a 
  +b2 a2 .  2  +b2 a2   −b2 a2  2  +b2 a2  
 
 2 2
b   −h2 b2   +h2 b2  
2 2
  −2h2 a2   −2h 
 −b h h2 2
a −2h 2 2
a −2h 2 2
a 
 
  −b a .   −2b a   +b a   −2b a  
2 2 2 2 2 2 2 2
−h2 a2 −2h2 b2  −2h2 b2  +h2 b2

  2 2     
ha −2h2 a2 h2 a2 h2 a2
 4  +b2 a2   +b2 a2   −2h2 b2  2  +b2 a2  
 
 2 2 2 2 2 2 2 2 
  +h2 b2  −2h2 b2   −2b 2a 2   −2h2 b 2  
 −2h a ha −2h a −2h a 
 
  +b a  4  +b a   −2b a  2  +b a  
2 2 2 2 2 2 2 2
 
 −2h2 b2   +h2 b2  +h2 b2  +h2 b2  

[Adi ] =     
h 2 2
a −2h 2 2
a h 2 2
a h 2 2
a 
 
  −2h b   −2b a  4  +b a  2  −2b a  
2 2 2 2 2 2 2 2
 
 −2b 2 2
a +h 2 2
b +h 2 2
b +h 2 2
b 
        2 2  
 ha2 2
−2h a 2 2
ha2 2
ha 
 
 2  +b2 a2  2  +b2 a2  2  −2b2 a2  4  +b2 a2  
−2h2 b2 +h2 b2 +h2 b2 +h2 b2
long
si b = h = a = n
tenemos:
117
 
4 0 0 0 −1 −1 −1 −1
 0 4 −1 −1 0 −1 0 −1 
 
 0 −1 4 −1 0 0 −1 −1 
 
klong 
 0 −1 −1 4 −1 0 0 −1 

AEk =
12n  −1 0 0 −1 4 −1 −1 0 

 −1 −1 0 0 −1 4 −1 0 
 
 −1 0 −1 0 −1 −1 4 0 
−1 −1 −1 −1 0 0 0 4

flEk = f (t) ŵl (x, y, z) dEk , l = 1, · · · , 8
Ek
a h b a h b
1
f1Ek = f (t) ŵ1 (x, y, z) dx dy dz = f (t) (b − x) (h − y) (a − z) dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= ((b − x) (h − y) (a − z)) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
a h b a h b
1
f2Ek = f (t) ŵ2 (x, y, z) dx dy dz = f (t) x (h − y) (a − z) dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= (x (h − y) z) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
a h b a h b
1
f3Ek = f (t) ŵ3 (x, y, z) dx dy dz = f (t) (b − x) y (a − z) dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= ((b − x) y (a − z)) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
a h b a h b
1
f4Ek = f (t) ŵ4 (x, y, z) dx dy dz = f (t) (b − x) (h − y) z dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= ((b − x) (h − y) z) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
a h b a h b
1
f5Ek = f (t) ŵ5 (x, y, z) dx dy dz = f (t) xy (a − z) dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= (xy (a − z)) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
118
a h b a h b
1
f6Ek = f (t) ŵ6 (x, y, z) dx dy dz = f (t) (b − x) yz dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= (xyz) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
a h b a h b
1
f7Ek = f (t) ŵ7 (x, y, z) dx dy dz = f (t) x (h − y) z dx dy dz
bha
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= (x (h − y) z) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
a h b a h b
xyz
f8Ek = f (t) ŵ8 (x, y, z) dx dy dz = f (t) bha
dx dy dz
0 0 0 0 0 0
a h b
f (t) f (t) 1 2 2 2 f (t)
= (xyz) dx dy dz = bha = (bha)
bha bha 8 8
0 0 0
f (t) (bha) T
f Ek = 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 .
8
long
si b = h = a = n
tenemos:
f (t) (long)3 T
f Ek = 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 .
8n3

qE
l
k ∩Γ2h
= q̂ (t) ( γ ŵl (x, y, z))|Ek ∩Γ2h d (Ek ∩ Γ2h ) , l = 1, · · · , 8.
Ek ∩Γ2h
Ek ∩ Γ2h tiene 6 posibles casos:

b1 , b2 , b4 , b7 , donde y = 0.
b3 , b5 , b6 , b8 , donde y = h.
b1 , b2 , b3 , b5 , donde z = 0.
b4 , b6 , b7 , b8 , donde z = a.
b1 , b3 , b4 , b6 , donde x = 0.
b2 , b5 , b7 , b8 , donde x = b.
- si b1 , b2 , b4 , b7 ∈ Ek ∩ Γ2h tenemos:
119
   

 a b 
 
 a b 


 
 
 (b−x)(h−y)(a−z) 

 q̂ (t) (γ ŵ1 (x, y, z))|y=0 dx dz 
 
 q̂ (t) bha dx dz 



 
 
 y=0 


 0 0 
 
 0 0 


 
 
 


 a b 
 
 a b 


 
 
 x(h−y)(a−z) 


 q̂ (t) (γ ŵ2 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=0 


 0 0 
 
 0 0 


 
 
 


 a b 
 
 a b 


 
 
 (b−x)y(a−z) 


 q̂ (t) (γ ŵ3 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=0 


 
 
 


 0 0 
 
 0 0 


 a b 
 
 a b 


 
 
 (b−x)(h−y)z



 q̂ (t) (γ ŵ4 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=0 

0 0 0 0
qEk ∩Γ2h = =

 a b 
 
 a b 


 
 
 xy(a−z) 


 q̂ (t) (γ ŵ5 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 

bha
y=0 


 0 0 
 
 0 0 


 
 
 


 a b 
 
 a b 


 
 
 (b−x)yz 


 q̂ (t) (γ ŵ6 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=0 


 
 
 



0 0 
 

0 0 


 a b 
 
 a b 


 
 
 x(h−y)z 


 q̂ (t) (γ ŵ7 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=0 


 
 
 


 0 0 
 
 0 0 


 a b 
 
 a b 


 
 
 xyz 


 q̂ (t) (γ ŵ8 (x, y, z))|y=0 dx dz 
 
 q̂ (t) dx dz 


 
 
 bha y=0 

0 0 0 0
 
 a b


 q̂(t) 


 (b − x) (a − z) dx dz 


 ba 


 



0 0 


 a b 


 
  1 2 2   


q̂(t)
x (a − z) dx dz 
  ba   1 

 ba 
 
 4
1 2 2 
 
 


 
 
 b a 
 
 1 


 0 0 
 
 4 
 
 


 0 
 
 0 
 
 0 


 
 
 1 2 2

 
 

a b q̂ (t) 4
b a q̂ (t) ba 1
= q̂(t) = =
 ba
 (b − x) z dx dz 
 ba   0 
 4   0 


 
 
 
 
 


 0 0 
 
 0 
 
 0 


 
 
 1 2 2   
 


 0 
 
 4
b a 
 
 1 


 
    

 0 
 0 0

 


 a b 


 q̂(t) 

 ba
 xz dx dz 


 


 

 0 0 
0
120
   

 a b 
 
 a b 


 
 
 (b−x)(h−y)(a−z) 

 q̂ (t) (γ ŵ1 (x, y, z))|y=h dx dz 
 
 q̂ (t) bha dx dz 



 
 
 y=h 


 0 0 
 
 0 0 


 
 
 


 a b 
 
 a b 


 
 
 x(h−y)(a−z) 


 q̂ (t) (γ ŵ2 (x, y, z))|y=h dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=h 


 0 0 
 
 0 0 


 
 
 


 a b 
 
 a b 


 
 
 (b−x)y(a−z) 


 q̂ (t) (γ ŵ3 (x, y, z))|y=hh dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=h 


 
 
 


 0 0 
 
 0 0 


 a b 
 
 a b 


 
 
 (b−x)(h−y)z



 q̂ (t) (γ ŵ4 (x, y, z))|y=h dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=h 

0 0 0 0
qEk ∩Γ2h = =

 a b 
 
 a b 


 
 
 xy(a−z) 


 q̂ (t) (γ ŵ5 (x, y, z))|y=h dx dz 
 
 q̂ (t) dx dz 


 
 

bha
y=h 


 0 0 
 
 0 0 


 
 
 


 a b 
 
 a b 


 
 
 (b−x)yz 


 q̂ (t) (γ ŵ6 (x, y, z))|y=h dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=h 


 
 
 



0 0 
 

0 0 


 a b 
 
 a b 


 
 
 x(h−y)z 


 q̂ (t) (γ ŵ7 (x, y, z))|y=h dx dz 
 
 q̂ (t) dx dz 


 
 
 bha
y=h 


 
 
 


 0 0 
 
 0 0 


 a b 
 
 a b 


 
 
 xyz 


 q̂ (t) (γ ŵ8 (x, y, z))|y=h dx dz 
 
 q̂ (t) dx dz 


 
 
 bha y=h 

0 0 0 0
 
 0 

 


 0 


 


 a b 


 q̂(t) 


 (b − x) (a − z) dx dz 


 ba 
    

 
  0   0 

 0 0 
 
 
 
 


 0 
 
 0 
 
 0 


 
 
 
 
 


 a b 
 
 1 2 2
ba 
 
 1 


 q̂(t) 
 
 4 
  
 
x (a − z) dx dz 0 0
= ba = q̂(t) q̂(t)ba
= 4

 0 0 

ba  14 b2 a2
 
  1 
 

 a b 
 
 1 2 2 
 
 


 
 
 4
ba 
 
 1 


 q̂(t) 
 
 
 
 


 (b − x) z dx dz 
 
 0 
 
 0 
 ba
 
  1 2 2    

 0 0 
 4
ba 1

 


 0 


 


 a b 


 q̂(t) 


 xz dx dz 

 ba 
0 0
121
   

 h b 
 
 h b 


 
 
 (b−x)(h−y)(a−z) 

 q̂ (t) (γ ŵ1 (x, y, z))|z=0 dx dy 
 
 q̂ (t) bha dx dy 



 
 
 z=0 


 0 0 
 
 0 0 


 
 
 


 h b 
 
 h b 


 
 
 x(h−y)(a−z) 


 q̂ (t) (γ ŵ2 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=0 


 0 0 
 
 0 0 


 
 
 


 h b 
 
 h b 


 
 
 (b−x)y(a−z) 


 q̂ (t) (γ ŵ3 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=0 


 
 
 


 0 0 
 
 0 0 


 h b 
 
 h b 


 
 
 (b−x)(h−y)z



 q̂ (t) (γ ŵ4 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=0 

0 0 0 0
qEk ∩Γ2h = =

 h b 
 
 h b 


 
 
 xy(a−z) 


 q̂ (t) (γ ŵ5 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 

bha
z=0 


 0 0 
 
 0 0 


 
 
 


 h b 
 
 h b 


 
 
 (b−x)yz 


 q̂ (t) (γ ŵ6 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=0 


 
 
 



0 0 
 

0 0 


 h b 
 
 h b 


 
 
 x(h−y)z 


 q̂ (t) (γ ŵ7 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=0 


 
 
 


 0 0 
 
 0 0 


 h b 
 
 h b 


 
 
 xyz 


 q̂ (t) (γ ŵ8 (x, y, z))|z=0 dx dy 
 
 q̂ (t) dx dy 


 
 
 bha z=0 

0 0 0 0
 
 h b 

 q̂(t) 


 (b − x) (h − y) dx dy 


 bh 


 



0 0 



h b 


 
    


q̂(t)
x (h − y) dx dy 
  1 2 2
bh   1 

 bh 
 
 4 
 
 


 
 
 1 2 2
bh 
 
 1 


 0 0 
 
 4 
 
 



h b 
 
 1 2 2
bh 
 
 1 


 q̂(t)

 
 4 
 
 

(b − x) y dx dy q̂ (t) 0 q̂ (t) bh 0
= bh = 1 2 2 =

 0 0 
 bh   4
bh  4   1 


 
 
 
 
 


 0 
 
 0   
 0 


 
 
 
 
 


 h b 
 
 0   
 0 

 q̂(t) 
     

 xy dx dy 
 0 0

 bh 


 



0 0 


 0 


 


 0 

 
0
- si b4 , b6 , b7 , b8 ∈ Ek ∩ Γ2h
122
   

 h b 
 
 h b 


 
 
 (b−x)(h−y)(a−z) 

 q̂ (t) (γ ŵ1 (x, y, z))|z=a dx dy 
 
 q̂ (t) bha dx dy 



 
 
 z=a 


 0 0 
 
 0 0 


 
 
 


 h b 
 
 h b 


 
 
 x(h−y)(a−z) 


 q̂ (t) (γ ŵ2 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=a 


 0 0 
 
 0 0 


 
 
 


 h b 
 
 h b 


 
 
 (b−x)y(a−z) 


 q̂ (t) (γ ŵ3 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=a 


 
 
 


 0 0 
 
 0 0 


 h b 
 
 h b 


 
 
 (b−x)(h−y)z



 q̂ (t) (γ ŵ4 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=a 

0 0 0 0
qEk ∩Γ2h = =

 h b 
 
 h b 


 
 
 xy(a−z) 


 q̂ (t) (γ ŵ5 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 

bha
z=a 


 0 0 
 
 0 0 


 
 
 


 h b 
 
 h b 


 
 
 (b−x)yz 


 q̂ (t) (γ ŵ6 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=a 


 
 
 



0 0 
 

0 0 


 h b 
 
 h b 


 
 
 x(h−y)z 


 q̂ (t) (γ ŵ7 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 
 bha
z=a 


 
 
 


 0 0 
 
 0 0 


 h b 
 
 h b 


 
 
 xyz 


 q̂ (t) (γ ŵ8 (x, y, z))|z=a dx dy 
 
 q̂ (t) dx dy 


 
 
 bha z=a 

0 0 0 0
 
 0 

 


 0 


 


 0 


 



h b 


 q̂(t) 
    

 (b − x) (h − y) dx dy 
  0   0 

 bh 
 
 
 
 


 
 
 0 
 
 0 


 0 0 
 
 
 
 


 0 
 
 0 
 
 0 


 
    
h b q̂ (t)  1 2 2
4
bh
 q̂ (t) bh 
1

= q̂(t) = =

 bh
(b − x) y dx dy 
 bh   0 
 4   0 


 
 
 1 2 2 
 
 


 0 0 
 
 4
bh 
 
 1 


 h b 
 
 1 2 2 
 
 


 
 
 bh 
 
 1 

 q̂(t) 
  4
1 2 2    

 x (h − y) dx dy 
 4
bh 1


bh 


 0 0 


 


 h b 


 q̂(t) 


 xy dx dy 

 bh 
0 0
- si b1 , b3 , b4 , b6 ∈ Ek ∩ Γ2h
123
   

 a h 
 
 a h 


 
 
 (b−x)(h−y)(a−z) 

 q̂ (t) (γ ŵ1 (x, y, z))|x=0 dy dz 
 
 q̂ (t) bha dy dz 



 
 
 x=0 


 0 0 
 
 0 0 


 
 
 


 a h 
 
 a h 


 
 
 x(h−y)(a−z) 


 q̂ (t) (γ ŵ2 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=0 


 0 0 
 
 0 0 


 
 
 


 a h 
 
 a h 


 
 
 (b−x)y(a−z) 


 q̂ (t) (γ ŵ3 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=0 


 
 
 


 0 0 
 
 0 0 


 a h 
 
 a h 


 
 
 (b−x)(h−y)z



 q̂ (t) (γ ŵ4 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=0 

0 0 0 0
qEk ∩Γ2h = =

 a h 
 
 a h 


 
 
 xy(a−z) 


 q̂ (t) (γ ŵ5 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 

bha
x=0 


 0 0 
 
 0 0 


 
 
 


 a h 
 
 a h 


 
 
 (b−x)yz 


 q̂ (t) (γ ŵ6 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=0 


 
 
 



0 0 
 

0 0 


 a h 
 
 a h 


 
 
 x(h−y)z 


 q̂ (t) (γ ŵ7 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=0 


 
 
 


 0 0 
 
 0 0 


 a h 
 
 a h 


 
 
 xyz 


 q̂ (t) (γ ŵ8 (x, y, z))|x=0 dy dz 
 
 q̂ (t) dy dz 


 
 
 bha x=0 

0 0 0 0
 
 a h


 q̂(t) 


 (h − y) (a − z) dy dz 


 ha 


 



0 0 


 0 


 
  1 2 2   


a h 
  ha   1 

 
 
 4 
 
 



q̂(t)
y (a − z) dy dz 
 
 0 
 
 0 


 ha 
 
 
 
 


 
 
 1 2 2
h a 
 
 1 


 0 0 
 
 4
1

 
 


a h q̂ (t) 4
h a2
2 q̂ (t) ha 1
= q̂(t) = =
 ha
 (h − y) z dy dz 
 ha   0 
 4  0 
 

 
 
 1 2 2 
 
 


 0 0 
 
 4
ha 
 
 1 


 
 
 
 
 


 0 
 
 0 
 
 0 

 
     

 a h 
 0 0

 q̂(t) 


 yz dy dz 


 ha 


 



0 0 


 0 

 
0
- si b2 , b5 , b7 , b8 ∈ Ek ∩ Γ2h
124
   

 a h 
 
 a h 


 
 
 (b−x)(h−y)(a−z) 

 q̂ (t) (γ ŵ1 (x, y, z))|x=b dy dz 
 
 q̂ (t) bha dy dz 



 
 
 x=b 


 0 0 
 
 0 0 


 
 
 


 a h 
 
 a h 


 
 
 x(h−y)(a−z) 


 q̂ (t) (γ ŵ2 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=b 


 0 0 
 
 0 0 


 
 
 


 a h 
 
 a h 


 
 
 (b−x)y(a−z) 


 q̂ (t) (γ ŵ3 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=b 


 
 
 


 0 0 
 
 0 0 


 a h 
 
 a h 


 
 
 (b−x)(h−y)z



 q̂ (t) (γ ŵ4 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=b 

0 0 0 0
qEk ∩Γ2h = =

 a h 
 
 a h 


 
 
 xy(a−z) 


 q̂ (t) (γ ŵ5 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 

bha
x=b 


 0 0 
 
 0 0 


 
 
 


 a h 
 
 a h 


 
 
 (b−x)yz 


 q̂ (t) (γ ŵ6 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=b 


 
 
 



0 0 
 

0 0 


 a h 
 
 a h 


 
 
 x(h−y)z 


 q̂ (t) (γ ŵ7 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 
 bha
x=b 


 
 
 


 0 0 
 
 0 0 


 a h 
 
 a h 


 
 
 xyz 


 q̂ (t) (γ ŵ8 (x, y, z))|x=b dy dz 
 
 q̂ (t) dy dz 


 
 
 bha x=b 

0 0 0 0
 
 0 

 


 a h 


 q̂(t) 


 (h − y) (a − z) dy dz 


 ha 


 


 0 0 
    

 0 
  0   0 

 
 
 
 
 


 0 
 
 1 2 2
ha 
 
 1 


 
 
 4 
 
 


 a h 
 
 0 
 
 0 


 q̂(t) 
 q̂ (t) 
 
 
 

y (a − z) dy dz 0 q̂ (t) ha 0
= ha = 1 2 2 =

 0 0 
 ha   4
h a 
 4  1 
 

 
 
 
 
 


 0 
 
 0 
 
 0 


 a h 
 
 1 2 2 
 
 


 
 
 ha 
 
 1 
 q̂(t)
 
  41 2 2    

 (h − y) z dy dz 
 4
ha 1


ha 


 0 0 


 


 a h 


 q̂(t) 


 yz dy dz 

 ha 
0 0
125
q̂ (t) ba T
qEk ∩Γ2h = 1 1 0 1 0 0 1 0 en b1 , b2 , b4 , b7 .
4
q̂ (t) ba T
qEk ∩Γ2h = 0 0 1 0 1 1 0 1 en b3 , b5 , b6 , b8 .
4
q̂ (t) bh T
qEk ∩Γ2h = 1 1 1 0 1 0 0 0 en b1 , b2 , b3 , b5 .
4
q̂ (t) bh T
qEk ∩Γ2h = 0 0 0 1 0 1 1 1 en b4 , b6 , b7 , b8 .
4
q̂ (t) ha T
qEk ∩Γ2h = 1 0 1 1 0 1 0 0 en b1 , b3 , b4 , b6 .
4
q̂ (t) ha T
qEk ∩Γ2h = 0 1 0 0 1 0 1 1 en b2 , b5 , b7 , b8 .
4
long
si b = h = a = n
tenemos:
q̂ (t) (long)2 T
q Ek ∩Γ2h
= 1 1 0 1 0 0 1 0 en b1 , b2 , b4 , b7 .
4n2 2
q̂ (t) (long) T
qEk ∩Γ2h = 0 0 1 0 1 1 0 1 en b3 , b5 , b6 , b8 .
4n2
q̂ (t) (long)2 T
qEk ∩Γ2h = 1 1 1 0 1 0 0 0 en b1 , b2 , b3 , b5 .
4n2
q̂ (t) (long)2 T
qEk ∩Γ2h = 0 0 0 1 0 1 1 1 en b4 , b6 , b7 , b8 .
4n2 2
q̂ (t) (long) T
qEk ∩Γ2h = 1 0 1 1 0 1 0 0 en b1 , b3 , b4 , b6 .
4n2 2
q̂ (t) (long) T
qEk ∩Γ2h = 0 1 0 0 1 0 1 1 en b2 , b5 , b7 , b8 .
4n2
2.10. Caso estacionario unidimensional.

Considerando el modelo discreto [2.38 − 2.40] expresamos el modelo discreto estacionario como:
Encuentre υ ∈ Kmh :
0 ≥ −Aυ· {η − υ} + f· {η − υ} , ∀η ∈ Kmh . (2.61)
Kmh = η ∈ Rmh | η > 0 en Ωh , γη|Γ1h = û sobre Γ1h . (2.62)

En 2.62 las restricciones deben ser introducidas al modelo 2.61, en el caso de la restricción η > 0 en
Ωh se cumple al seleccionar los valores de la fuente interna, (ver simulaciones numéricas,
figuras 2.11-
υ η
2.14), y para la restricción γη|Γ1h = û sobre Γ1h se requiere seleccionar υ = yη=
û û
substituyendo en 2.61,

υ η υ η υ
0 ≥ −A · − + f· −
û û û û û
126
$ %
υ 0 η−υ η−υ
0 ≥ −A + · + f·
0 û 0 0

υ η−υ 0 η−υ η−υ
0 ≥ −A · −A · + f·
0 0 û 0 0

υ η−υ 0 η−υ η−υ
0 ≥ −A · −A · + f·
0 0 û 0 0
$ %
υ η−υ 0 η−υ
0 ≥ −A · + −A +f ·
0 0 û 0
$ %
0
0 ≥ −Aυ· {η − υ} + −A + f · {η − υ}
û
* $ %+
0
0 ≥ −Aυ + −A +f · {η − υ}
û
$ %
0
−Aυ + −A +f =0
û
$ %
0
aquí Aef = A es la matriz de conductividad efectiva, y fef = −A +f es el vector efectivo
û
de cargas. El problema estacionario resulta ser:
Encuentre υ ∈ Rm ,h
:
Aef υ = fef . (2.63)
Iniciamos obteniendo la solución del caso estacionario, de acuerdo a diversas discretizaciones es-
paciales y también diferentes valores de fuente interna, lo anterior para tener como referencia las
soluciones estacionarias, ya que el problema evolutivo revisado posteriormente es independiente de
la condición inicial tiene como límite su estado estacionario.
Además, por los algoritmos evolutivos requerimos tener el valor del parámetro w, el cual puede ser
calculado teniendo las soluciones estacionarias correspondientes.
Consideremos
W las discretizaciones: n [elementos] = 5, 10, 15, · · · , 100 y la fuente interna con valores:
f m3 = -150000, -125000, -100000, -50000, -10000, -5000, -1000, -500, 0, 500, 1000, 5000, 10000,
50000, 100000, 125000, 150000.
En lo anterior se tienen que f m W
3 ≥ −150000 hace que la suposición u > 0 sea satisfecha.
Las condiciones sobre la frontera son: son: û(0) = 300 [K] y û(1.0) = 500 [K] .
Se resuelven los casos correspondientes a las 20 discretizaciones y 17 valores de fuente interna, y se
presentan los casos n= 5, 35, 65 y 100 los cuales gráficamente son:
127
Fig. 2.11. Soluciones Estacionarias , para n=5.
128
f
Calculamos el parámetro ω = para cada solución y obtenemos la gráfica de comportamiento
A u
del parámetro ω :
f
Fig. 2.15. Parámetro ω =
A u
Seleccionamos un caso de fuente interna f = 10000, con el cual se harán las simulaciones de los
algoritmos de integración temporal, en este caso tenemos los valores del parámetro ω para cada
discretización n son dados en la tabla:
n ω n ω n ω n ω
5 0.194402686 30 0.06450989 55 0.047180411 80 0.038991624
10 0.118720884 35 0.059531964 60 0.045131338 85 0.037811924
15 0.093700065 40 0.055557625 65 0.043328413 90 0.036733300
20 0.079998357 45 0.052288477 70 0.041725932 95 0.035742086
25 0.071004728 50 0.049537442 75 0.040289264 100 0.034827072.
129
2.10.1. Código computacional.
Escribimos el código computacional en lenguaje Fortran77.
C23456789012345678901234567890123456789012345678901234567890123456789012
CJorge Sanchez Macias febrero 2007.
C programa para solucion estacionaria y dobleus w
C Caso: ’’unidimensional’’
C
PROGRAM estacionario
Implicit none
Integer n,i,j
Double Precision cc,den,cevol,are,lon,fint,
* kk(101,101),ff(101),kloc(2,2),efe(2),
* kkk(99,99),uupaso1(101),ddd(99),beta,unorm,fnorm,
* fff(99),kkmod(99,2),cf(2),otro(101),dobleus
character*12 salida
C
C n asigna el numero de elementos
do n=5,100,5
C
C asigna nombre a archivo de salida.
select case (n)
case (5)
salida=’’est051d.dat’’
case (10)
case (15)
case (20)
case (25)
case (30)
case (35)
case (40)
case (45)
case (50)
130
case (55)
case (60)
case (65)
case (70)
case (75)
case (80)
case (85)
case (90)
case (95)
case (100)
case default
salida=’’default.dat’’
end select
write(*,10) salida
C datos: calor especifico a volumen constante, densidad,
C conductividad termica, area, longitud.fuente interna
cevol= 480.0d0
den=7850.0d0
cc= 50.0d0
are=0.001d0
lon= 1.0d0
C
fint= 10000.0d0
C
C definicion matrices locales
kloc(1,1)= 1.0d0*cc*are*n/lon
kloc(1,2)= -1.0d0*cc*are*n/lon
efe(1)= 1.0d0*fint*are*lon/2.0d0/n
131
C limpia matrices globales
do i=1,101,1
do j=1,101,1
kk(i,j)= 0.0d0
enddo
ff(i)= 0.0d0
otro(i)= 0.0d0
enddo
C limpia matrices reducidas
do i=1,99,1
do j=1,99,1
kkk(i,j)= 0.0d0
enddo
fff(i)= 0.0d0
kkmod(i,1)= 0.0d0
kkmod(i,2)= 0.0d0
enddo
C ensamble de matrices y vectores globales
do i=1,n,1
kk(i,i)= kk(i,i) + kloc(1,1)
kk(i,i+1)= kk(i,i+1) + kloc(1,2)
kk(i+1,i)= kk(i+1,i) + kloc(2,1)
kk(i+1,i+1)= kk(i+1,i+1) + kloc(2,2)
ff(i)= ff(i) + efe(1)
ff(i+1)= ff(i+1) + efe(2)
enddo
C matrices reducidas
do i=1,(n-1),1
do j=1,(n-1),1
kkk(i,j)= kk(i+1,j+1)
enddo
fff(i)= ff(i+1)
kkmod(i,1)= kk(i+1,1)
kkmod(i,2)= kk(i+1,n+1)
enddo
C normas de matrices
do i=1,99,1
ddd(i)= 0.0d0
enddo
call jacobi(kkk,(n-1),99,ddd)
call eigsrt(ddd,(n-1),99)
beta= ddd(1)
132
do i=1,(n-1),1
do j=1,(n-1),1
enddo
enddo
C condiciones de frontera
cf(1)=300
cf(2)=500
C vector efectivo de cargas
do i=1,(n-1),1
fff(i)= fff(i)-kkmod(i,1)*cf(1)-kkmod(i,2)*cf(2)
enddo
C calcula norma vector efectivo de cargas
call enor(99,(n-1),fff,fnorm)
C actualiza matriz de conductividad A
do i=1,99,1
do j=1,99,1
kk(i,j)= kkk(i,j)
enddo
ff(i)= fff(i)
enddo
C descomposicion cholesky
call choldc(kk,(n-1),101,otro)
C resuelve
call subpaso(101,(n-1),kk,ff,uupaso1)
C calcula norma vector solucion
call enor(101,(n-1),uupaso1,unorm)
dobleus=fnorm/beta/unorm
C guarda solucion
open(3, FILE = salida)
write(3,*)cf(1)
do i=1,(n-1),1
write(3,*)uupaso1(i)
enddo
write(3,*)cf(2)
write(3,*)dobleus
write(3,*)fnorm
write(3,*)beta
write(3,*)unorm
write(3,*)n
write(3,*)fint
close(3)
133
enddo
C
10 format(A15)
15 format(I3)
20 format(6F12.6)
25 format(4F12.6)
30 format(4F12.6) C
END
C
SUBROUTINE jacobi(a,n,np,d)
C version modificada de programa original de numerical recipes fortran
INTEGER n,np,nrot,NMAX
Double Precision a(np,np),d(np),v(np,np)
PARAMETER (NMAX=500)
INTEGER i,ip,iq,j
Double Precision c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
do 12 ip=1,n
do 11 iq=1,n
v(ip,iq)=0.0d0
11 continue
v(ip,ip)=1.0d0
12 continue
do 13 ip=1,n
b(ip)=a(ip,ip)
d(ip)=b(ip)
z(ip)=0.0d0
13 continue
nrot=0
do 24 i=1,50
sm=0.
do 15 ip=1,n-1
do 14 iq=ip+1,n
sm=sm+dabs(a(ip,iq))
14 continue
15 continue
if(sm.eq.0.)return
if(i.lt.4)then
tresh=0.2d0*sm/n**2
else
tresh=0.
endif
do 22 ip=1,n-1
134
do 21 iq=ip+1,n
g=100.*dabs(a(ip,iq))
if((i.gt.4).and.(dabs(d(ip))+
*g.eq.dabs(d(ip))).and.(dabs(d(iq))+g.eq.dabs(d(iq))))then
a(ip,iq)=0.0d0
else if(dabs(a(ip,iq)).gt.tresh)then
h=d(iq)-d(ip)
if(dabs(h)+g.eq.dabs(h))then
t=a(ip,iq)/h
else
theta=0.5d0*h/a(ip,iq)
t=1.0d0/(dabs(theta)+dsqrt(1.0d0+theta**2))
if(theta.lt.0.)t=-t
endif
c=1.0d0/sqrt(1+t**2)
s=t*c
tau=s/(1.0d0+c)
h=t*a(ip,iq)
z(ip)=z(ip)-h
z(iq)=z(iq)+h
d(ip)=d(ip)-h
d(iq)=d(iq)+h
a(ip,iq)=0.0d0
do 16 j=1,ip-1
g=a(j,ip)
h=a(j,iq)
a(j,ip)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
16 continue
do 17 j=ip+1,iq-1
g=a(ip,j)
h=a(j,iq)
a(ip,j)=g-s*(h+g*tau)
17 continue
do 18 j=iq+1,n
g=a(ip,j)
h=a(iq,j)
a(iq,j)=h+s*(g-h*tau)
18 continue
do 19 j=1,n
135
g=v(j,ip)
h=v(j,iq)
v(j,ip)=g-s*(h+g*tau)
v(j,iq)=h+s*(g-h*tau)
19 continue
nrot=nrot+1
endif
21 continue
22 continue
do 23 ip=1,n
b(ip)=b(ip)+z(ip)
d(ip)=b(ip)
z(ip)=0.0d0
23 continue
24 continue
pause ’too many iterations in jacobi’
return
END
C
SUBROUTINE eigsrt(d,n,np)
INTEGER n,np
Double Precision d(np),v(np,np)
INTEGER i,j,k
Double Precision p
do 13 i=1,n-1
k=i
p=d(i)
do 11 j=i+1,n
if(d(j).ge.p)then
k=j
p=d(j)
endif
11 continue
if(k.ne.i)then
d(k)=d(i)
d(i)=p
do 12 j=1,n
p=v(j,i)
v(j,i)=v(j,k)
v(j,k)=p
12 continue
136
endif
13 continue
return
END
C
SUBROUTINE subpaso(np,n,m1,vect,sact)
Integer n,np,i,k
Double Precision m1(np,np),vect(np),sact(np),c(n)
C
C resuelve c por triangular inferior: m1*c=b
do i=1,n,1
c(i)=vect(i)
do k=1,i-1,1
c(i)=c(i)-m1(i,k)*c(k)
enddo
c(i)=c(i)/m1(i,i)
enddo
C resuelve sact por triangular superior: m1*sact=c
do i=1,np,1
sact(i)=0.0d0
enddo
do i=n,1,-1
sact(i)=c(i)
do k=i+1,n,1
sact(i)=sact(i)-m1(i,k)*sact(k)
enddo
sact(i)=sact(i)/m1(i,i)
enddo
return
END
C
SUBROUTINE choldc(a,n,np,p)
INTEGER n,np
Double Precision a(np,np),p(np)
INTEGER i,j,k
Double Precision sum
do 13 i=1,n
do 12 j=i,n
sum=a(i,j)
do 11 k=i-1,1,-1
sum=sum-a(i,k)*a(j,k)
137
11 continue
if(i.eq.j)then
p(i)=dsqrt(sum)
else
a(j,i)=sum/p(i)
endif
12 continue
13 continue
do 15 i=1,n
do 14 j=i+1,n
a(i,j)=a(j,i)
14 continue
a(i,i)=p(i)
15 continue
return
END
C
SUBROUTINE enor(np,n,um1,er)
Double Precision um1(np),er,su1
Integer i
su1= 0.0d0
do i=1,n,1
su1= su1 + (um1(i)*um1(i))
enddo
er=dsqrt(su1)
END
C
TAREA 6. Ejecute el programa anterior para los siguientes casos:

n = (15 + n.l.) f = 22000
n = (40 + n.l.) f = 22000
n = (65 + n.l.) f = 22000
los demás datos y condiciones de frontera se mantienen de acuerdo al ejemplo unidimensional.
Imprima el archivo de salida (solución y datos} de cada caso.
Genere la grafica de las tres soluciones juntas.
Analice la gráfica y discuta brevemente lo más relevante o importante de las soluciones (de media a
una cuartilla.).
2.11. Caso evolutivo unidimensional.

Considerando el modelo discreto [2.38 − 2.40] tenemos:
138
υ (0) = υ 0 en Ωh ,
encuentre υ ∈ Ktmh :
•
0 ≥ − M υ · {η − υ} − Aυ· {η − υ} + f · {η − υ} , ∀η ∈ Ktmh .
Ktmh = η : R → Rmh | η (t) > 0 en Ωh , γη (t)|Γ1h = û sobre Γ1h .
En el conjunto de temperaturas admisibles Ktmh tenemos restricciones que deben ser introducidas al
modelo, en el caso de la restricción η (t) > 0 en Ωh se cumple al seleccionar los valores de la fuente
interna, (ver simulaciones numéricas estacionarias), y para la restricción γη (t)|Γ1h = û sobre Γ1h se
-
•
υ η • υ , substituyendo,
requiere seleccionar υ = yη= , ademas tenemos que υ=
û û 0
-
•
υ η υ υ η υ η υ
0 ≥ −M · − −A · − + f· −
0 û û û û û û û
- $ %
•
υ η−υ υ 0 η−υ η−υ
0 ≥ −M · −A + · + f·
0 0 0 û 0 0
-
•
υ η−υ υ η−υ 0 η−υ η−υ
0 ≥ −M · −A · −A · + f·
0 0 0 0 û 0 0
-
•
υ η−υ υ η−υ 0 η−υ η−υ
0 ≥ −M · −A · −A · + f·
0 0 0 0 û 0 0
* - + $ %
•
υ υ η−υ 0 η−υ
0 ≥ −M −A · + −A +f ·
0 0 0 û 0
$ %
• 0
0 ≥ −Mυ − Aυ · {η − υ} + −A + f · {η − υ}
û
* $ %+
• 0
0 ≥ −Mυ − Aυ + −A +f · {η − υ}
û
Reescribimos el modelo discreto como:

υ (0) = υ 0 en Ωh ,
•

,h ,h
encuentre υ ∈ Ktm : 0 ≥ −Mef υ −Aef υ + fef · {η − υ} , ∀η ∈ Ktm .
,h
Ktm = ,h
η : R → Rm .
139
o equivalentemente tenemos:
υ (0) = υ 0 en Ωh , (2.64)
,h •
encuentre υ ∈ Ktm : 0 ∈ −Mef υ −αAef υ − (1 − α) Aef υ + fef (2.65)
aquí Mef = M es la matriz de acumulación efectiva, Aef = A es la matriz de conductividad efectiva,

$ %
0
y fef = −A + f es el vector efectivo de cargas, además se ha introducido la descomposicón
û
convexa del término Aef υ con α ∈ (0, 1) .
El modelo 2.64 − 2.65 es la forma inicial de diversos procesos de integración temporal numérica como
los que aplicaremos, definimos algunos valores relacionados con las simulaciones numéricas.
Las condiciones sobre la frontera utilizadas son: û(0) =
300 Ky û(1.0) = 500 K.

 400  
Se considera como condición inicial, υ ∈ R : υ =
0 n−1 0 .
.
 .
 400  
Para el error relativo normalizado ε son utilizados los valores: 10−2 , 10−3 , 10−4 , · · · , 10−12 .
El número de elementos finitos n en los casos simulados son: 5, 10, 15, 20, · · · , 100.
Se procede a realizar las simulaciones numéricas, tomando en cuenta cuatro casos para el término ω
f
con valores: 0.0, , 0.5 y 1.0.
A u
Los algoritmos numéricos tratados para la resolución de 2.64 − 2.65 son: Semi-implícito de Euler,
Douglas-Rachford, Peaceman-Rachford, Esquema-θ.
Para cada esquema se presentan los resultados del comportamiento numérico de los algoritmos de
acuerdo a los datos declarados anteriormente y el en caso del esquema semi-implícito de Euler se
presenta la deducción del algoritmo y la codificación correspondiente.
2.11.1. Semi − Implı́cito de Euler.
Proposición. Bajo la simplificación de notación del problema evolutivo de primer orden discreto-
espacial (2.64 − 2.65) , donde, el superíndice corresponde a la evaluación del campo en tiempo dis-
creto, esto es, um = u (tm ) , um+1 = u (tm+1 ) , y además para las matrices del problema tenemos
M = Mef , A = Aef y f = f ef , escribimos el algoritmo Semi − implı́cito de Euler como:
140


 0.Determine β y µ : β = A y µ = M .

  Seleccione ω tal que: 0 ≤ ω ≤ 1, o



 f

  Determine ω como: ω = , donde u es solución de Au = f.



 A u

 Dado un error relativo normalizado límite ε suficientemente pequeño fijo.



 Seleccione ℓ tal que: 0 < ε < ℓ < 2$(1 + ω) .

 %

 2ℓ µ
 Determine k tal que cumpla: k = .
2 (1 + ω) − ℓ β (2.66)

 1. Dado u0 .



 2. Conocido um , para m ≥ 0.



 3. Calcular um+1 :

 −1



 um+1 = M + 12 kA M − 12 kA um + kf m+1 .



 um+1 − um

 4. Si ≤ ε, entonces terminar,

 um

Si no, entonces hacer m = m + 1 e ir al paso 2.
Demostración.
El proceso de demostración se realiza en tres partes: la primera establece el algoritmo (pasos 1 a 3) ,
la segunda caracteriza el error relativo normalizado (paso 4) y la tercera explicita la parametrización
del paso de tiempo (paso 0) .
Primera Parte.
Para problema evolutivo de primer orden discreto-espacial (2.65) su discretización temporal de
acuerdo al esquema Semi-Implícito de Euler es:
um+1 − um
0∈M + αAum + (1 − α) Aum+1 − f m+1 , (2.67)
k
donde k = tm+1 − tm > 0, con tm+1 y tm tiempos discretos para el número del paso de tiempo (la
iteración), m + 1, con m ≥ 0.
Nota 1. Una variante del esquema Semi-Implícito de Euler (2.65) es:
um+1 − um
0∈M + αAum+1 + (1 − α) Aum − f m+1 . (2.68)
k
A partir de (2.67) procedemos a explicitar el algoritmo correspondiente:

0 ∈ Mum+1 − Mum + kαAum + k (1 − α) Aum+1 − kf m+1
0 ∈ Mum+1 + k (1 − α) Aum+1 − Mum + kαAum − kf m+1
0 ∈ (M + k (1 − α) A) um+1 − (M − kαA) um − kf m+1
obtenemos el algoritmo:

um+1 = (M + k (1 − α) A)−1 (M − kαA) um + kf m+1 . (2.69)
141
Nota 2. De la variante (2.68) , obtenemos equivalentemente:

um+1 = (M + kαA)−1 (M − k (1 − α) A) um + kf m+1 . (2.70)
En (2.69) los términos (M + k (1 − α) A)−1 y (M − kαA) son contracciones bajo las condiciones
adecuadas del parámetro k, Gabay [Ga 1983] y donde la máxima velocidad de convergencia se obten-
drá cuando los términos k (1 − α) y kα de los operadores (M + k (1 − α) A)−1 y (M − kαA) tengan
los máximos valores posibles.
Graficando, α y (1 − α) tenemos,
0.8
0.6
α, (1 − α)
0.4
0.2
0 0 0.2 0.4 0.6 0.8 1
α
F ig.1.16.
de donde observamos que los términos k (1 − α) y kα tendrán valores máximos si: k (1 − α) = kα
de donde obtenemos α = 12 , sustituyendo en (2.69) y simplificando:
$ %−1 $$ % %
m+1 1 1 m m+1
u = M + kA M − kA u + kf . (2.71)
2 2
−1
Por los operadores presentes M + 12 kA y M − 12 kA en 2.71 tenemos que es convergente,
tomemos: limm→∞ um = u, limm→∞ um+1 = u, limm→∞ f m+1 = f entonces
$ %−1 $$ % %
m+1 1 1 m m+1
lim u = lim M + kA M − kA u + kf
m→∞ m→∞ 2 2
$ %−1 $$ % %
1 1 m m+1
= M + kA M − kA lim u + k lim f
2 2 m→∞ m→∞
$ %−1 $$ % %
1 1
u = M + kA M − kA u + kf
2 2
$ % $ %
1 1
M + kA u = M − kA u + kf
2 2
$ % $ %
1 1
M + kA u− M − kA u = kf
2 2
142
Au = f (2.72)
por lo tanto la solución um+1 de 2.71 converge hacia u, la cual es solución del problema estacionario
2.72.
Nota 3. En el caso del algoritmo (2.70) , al hacer la caracterización de α tenemos el mismo resultado:
α = 12 y el algoritmo resultante coincide con (2.71) .
Segunda Parte.
Expresamos (2.71) en términos de operadores resolventes, procediendo a explicitar um+1 − um :
$$ % %
m+1
1
k 1 m m+1
u = JM,A
2
M − kA u + kf
2
1
k −1
donde JM,A
2
= M + 12 kA
$ %
m+1
1
k 1 1
k m+1
u = JM,A M − kA um + JM,A
2 2
kf
2
$ $ % %
1
k 1 1
k m+1
um+1 − um = JM,A M − kA − I um + JM,A
2 2
kf
2
$$ % −1 %
1
k 1 1
k 1
k m+1
um+1 − um = JM,A
2
M − kA + JM,A 2
(−I) um + JM,A 2
kf
2
$$ % $ %%
1
k 1 1 1
k m+1
um+1 − um = JM,A
2
M − kA − M + kA um + JM,A
2
kf
2 2
1
k
um+1 − um = JM,A
2
(−kA) um + kf m+1
um+1 − um
obtenemos la expresión , considerando um = 0 si m ≥ 0 :
um
. m+1 . . 1 . . .
.u m. . 2k .
−u ≤ .JM,A . −kA um + .kf m+1 .
. 1 . $ %
. 2k . m kf m+1
≤ .JM,A . u −kA +
um
um+1 − um . . $ f m+1
%
. 12 k .
≤ .JM,A . k A +
um um
f m+1 f f f m+1
considerando limm→∞ = , de 2.72 tenemos ≤ A entonces = ω A ,
um u u um
donde 0 ≤ ω ≤ 1, y concluimos
um+1 − um .. 12 k .
.
≤ . J M,A . k (1 + ω) A . (2.73)
um
143
Definamos dos escalares estrictamente positivos, ε y ℓ, tales que 0 < ε < ℓ, con los cuales escribimos
la expresión anterior como:


 um+1 − um
≤ ε,
. 1 . u
m
donde 0 < ε < ℓ. (2.74)

 . k .
.JM,A
2
. k (1 + ω) A = ℓ,
Seleccionamos como proceso de detención del algoritmo (2.71) el error relativo normalizado dado por
(2.74)1 y donde el error relativo normalizado límite ε suficientemente pequeño es dado y fijo.
Tercera Parte.
Definamos β = A > 0 y µ = M > 0, y procedemos analizar (2.74)2 :
. 1 .
. 2k .
.JM,A . k (1 + ω) A = ℓ
1 . 1 .
. 2k .
como 1 ≤ .JM,A . , tenemos:
µ + 2 kβ
$ %
1
k (1 + ω) β ≤ ℓ
µ + 12 kβ
$ %
1
ℓ µ + kβ − k (1 + ω) β ≥ 0
2
$ %
1
ℓ − (1 + ω) βk + ℓµ ≥ 0
2
y establecemos
φ (k) = a k + b ≥ 0, (2.75)

donde a = 12 ℓ − (1 + ω) β y b = ℓµ.
Como ℓ > 0 y µ > 0 ⇒ b = ℓµ > 0, si a < 0 entonces φ (k1 ) = 0 donde k1 = − ab > 0.
De a < 0 tenemos ℓ < 2 (1 + ω) y tomando el valor de k = k1 , tenemos que (2.74) toma la forma:


 um+1 − um
 ≤ ε,
$ u
m
% (2.76)

 2ℓ µ
 k= ,
2 (1 + ω) − ℓ β
f
donde, 0 ≤ ω ≤ 1 ó ω = y además.0 < ε < ℓ < 2 (1 + ω) , concluyendo así la parame-
A u
trización de k y por lo tanto se cumple
144


 0.Determine β y µ : β = A y µ = M .

  Seleccione ω tal que: 0 ≤ ω ≤ 1, o



 f

  Determine ω como: ω = , donde u es solución de Au = f .



 A u

 Dado un error relativo normalizado límite ε suficientemente pequeño fijo.



 Seleccione ℓ tal que: 0 < ε < ℓ < 2$(1 + ω) .

 %

 2ℓ µ
 Determine k tal que cumpla: k = .
2 (1 + ω) − ℓ β

 1. Dado u0 .



 2. Conocido um , para m ≥ 0.



 3. Calcular um+1 :

 −1



 um+1 = M + 12 kA M − 12 kA um + kf m+1 .



 um+1 − um

 4. Si ≤ ε, entonces terminar,

 um


Presentamos los resultados de diferentes simulaciones numéricas.
Fig. 2.17
145
Fig. 2.18
Fig. 2.19
146
Fig. 2.20
En las figuras 2.16 - 2.20 tenemos el comportamiento numérico de la triangulación del dominio
n = 5, .., 100, observamos para el error relativo normalizado ε = 1E − 8 del algoritmo semi-implícito
f
de Euler para los valores ensayados de ω = 0, .0.5, 1.0 y .
A u
Fig.2.21
147
En la figura 2.21 tenemos la simulación de la triangulación n = 100, ω = 0, ε = 1E −8, con parámetro
l = 1.974 equivalente a un tamaño de paso en el tiempo k = 141.1667seg.
Comparando las triangulaciones n=100, n=50 y n=5, tenemos la tabla de valores:
w=0 n = 100 n = 50 n=5

k [seg.] [horas] k [seg.] [horas] k [seg.] [horas]
ε iter iter iter
1E − 2 356 126.1734 12.4772 184 272.1033 13.9075 20 3097.1920 15.4098
1E − 3 479 133.2611 17.7311 245 272.1033 18.5181 26 3167.7600 20.0327
1E − 4 604 133.2611 22.3583 307 272.1033 23.2044 31 3097.1920 24.9515
1E − 5 732 141.1667 27.0964 370 272.1033 27.9662 36 3097.1920 29.6219
1E − 6 856 141.1667 33.5663 434 286.4623 32.8036 42 3167.7600 34.5589
1E − 7 980 141.1667 38.4287 496 286.4623 39.4681 47 3097.1920 39.5414
1E − 8 1098 141.1667 43.2911 559 286.4623 44.4812 52 3097.1920 43.7479
1E − 9 1229 141.1667 48.1927 622 286.4623 49.4943 58 3097.1920 48.7957
1E − 10 1354 141.1667 53.0944 686 286.4623 54.5870 63 3097.1920 53.0023
1E − 11 1480 141.1667 58.0352 749 286.4623 59.6001 68 3097.1920 58.5025
1E − 12 4339 50.4520 60.8087 855 286.4623 68.0348 74 3097.1920 63.6645
Observamos que para las tres triangulaciones comparadas tenemos valores de tiempo de llegada al
estado estacionario cercanos entre sí.
Código computacional. Tenemos el código en lenguaje Fortran77.
C23456789012345678901234567890123456789012345678901234567890123456789012
C Jorge Sanchez Macias Mayo 2007
C Programa final para esquema Semi-Implicito de Euler
Program nuu
Implicit none
Integer iz1,iz2,n,i,j,anter,cont
Double Precision cevol,den,cc,are,lon,fint,distan,kloc(2,2),
* mloc(2,2),efe(2),kk(101,101),mm(101,101),ff(101),kkk(99,99),
* mmm(99,99),fff(99),ddd(99),kkmod(99,2),cf(2),otro(101),beta,mu,
* ele,kas,estac(101),ern,eab,uuini(101),uufin(101),dobleus
Logical out1,otravez
character*12 salida
C
C datos: calor especifico a volumen constante, densidad,
C conductividad te}rmica, area, longitud.fuente interna
cevol= 480.0d0
den= 7850.0d0
cc= 50.0d0
are= 0.001d0
lon= 1.0d0
148
fint= 10000.0d0
open(13,FILE=’’sienuu4.res’’)
C datos: calor especifico a volumen constante,densidad,
C conductividad termica, area, longitud.
write(13,*)’’datos: calor especifico a volumen constante,’’
write(13,*)’’densidad,conductividad termica,area,longitud,f.int.’’
write(13,*)cevol
write(13,*)den
write(13,*)cc
write(13,*)are
write(13,*)lon
write(13,*)fint
C
do iz1=8,8,1
distan=(1d0/10d0)**iz1
write(*,*)’’distancia= ’’,distan
do iz2=400,400,50
write(*,*)’’condicion inicial= ’’,iz2
do n=5,100,5
write(*,*)n
C asigna nombre a archivo de solucion estacionaria.
select case (n)
case (5)
case (10)
case (15)
case (20)
case (25)
case (30)
case (35)
case (40)
case (45)
case (50)
149
case (55)
case (60)
case (65)
case (70)
case (75)
case (80)
case (85)
case (90)
case (95)
case (100)
case default
salida=’’default.dat’’
end select
C lee solucion estacionaria y w
open(8,FILE=salida)
do i=1,(n+1),1
read(8,*)estac(i)
enddo
read(8,*)dobleus
close(8)
do i=(n+2),101,1
estac(i)= 0.0d0
enddo
C en el caso w><|f|/(|A||u|) se asigna el valor aqui
dobleus= 1.0d0
C asignar limite superior de lsup= 2(1+w)
ele= 3.960d0
C definicion matrices locales
mloc(1,1)= 2.0d0*cevol*den*are*lon/6.0d0/n
150
C limpia matrices globales
do i=1,101,1
do j=1,101,1
kk(i,j)= 0.0d0
mm(i,j)= 0.0d0
enddo
ff(i)= 0.0d0
otro(i)= 0.0d0
enddo
C limpia matrices reducidas
do i=1,99,1
do j=1,99,1
kkk(i,j)= 0.0d0
mmm(i,j)= 0.0d0
enddo
fff(i)= 0.0d0
kkmod(i,1)= 0.0d0
kkmod(i,2)= 0.0d0
enddo
C ensamble de matrices y vectores globales
do i=1,n,1
kk(i,i)= kk(i,i) + kloc(1,1)
kk(i,i+1)= kk(i,i+1) + kloc(1,2)
kk(i+1,i)= kk(i+1,i) + kloc(2,1)
kk(i+1,i+1)= kk(i+1,i+1) + kloc(2,2)
mm(i,i)= mm(i,i) + mloc(1,1)
mm(i,i+1)= mm(i,i+1) + mloc(1,2)
mm(i+1,i)= mm(i+1,i) + mloc(2,1)
mm(i+1,i+1)= mm(i+1,i+1) + mloc(2,2)
ff(i)= ff(i) + efe(1)
ff(i+1)= ff(i+1) + efe(2)
enddo
C matrices reducidas
do i=1,(n-1),1
do j=1,(n-1),1
151
mmm(i,j)= mm(i+1,j+1)
enddo
fff(i)= ff(i+1)
kkmod(i,1)= kk(i+1,1)
kkmod(i,2)= kk(i+1,n+1)
enddo
C normas de matrices
do i=1,99,1
ddd(i)= 0.0d0
enddo
call jacobi(kkk,(n-1),99,ddd)
beta= ddd(1)
do i=1,99,1
ddd(i)= 0.0d0
enddo
call jacobi(mmm,(n-1),99,ddd)
mu= ddd(1)
do i=1,(n-1),1
do j=1,(n-1),1
mmm(i,j)= mm(i+1,j+1)
enddo
enddo
C condiciones de frontera
cf(1)=300d0
cf(2)=500d0
C vector efectivo de cargas
do i=1,(n-1),1
fff(i)= kkmod(i,1)*cf(1)+kkmod(i,2)*cf(2)-fff(i)
enddo
C ciclo para ele
out1=.true.
anter=1000000
do while(out1)
ele= ele - 0.01d0
kas=((2.0d0*ele)/((2.0d0*(1.0d0+dobleus))-ele))*(mu/beta)
C limpia matrices y vectores
do i=1,101,1
do j=1,101,1
kk(i,j)= 0.0d0
152
mm(i,j)= 0.0d0
enddo
ff(i)= 0.0d0
enddo
C construye matrices y vectores para el esquema semi-implicito de euler
do i=1,(n-1),1
do j=1,(n-1),1
kk(i,j)= mmm(i,j)+(kas/2.0d0)*kkk(i,j)
mm(i,j)= mmm(i,j)-(kas/2.0d0)*kkk(i,j)
enddo
ff(i)= -(kas)*fff(i)
enddo
C descomposicion cholesky
call choldc(kk,(n-1),101,otro)
C asigna valor inicial
do i=1,(n-1),1
uuini(i)= iz2*1.0d0
enddo
do i=n,101,1
uuini(i)= 0.0d0
enddo
otravez=.TRUE.
ern= 1.0d0
cont= 0
do while (otravez)
cont= cont + 1
C resuelve un subpaso
call subpaso(101,(n-1),kk,mm,uuini,ff,uufin)
C calcula distancia entre soluciones consecutivas
call ernor(101,(n-1),uufin,uuini,ern)
C actualiza solucion
do i=1,101,1
uuini(i)= uufin(i)
enddo
C compara distancia entre soluciones consecutivas
if (ern.LE.distan)then
otravez=.FALSE.
endif
enddo
C calcula error absoluto
eab=0.0d0
do i=1,(n-1),1
153
eab= eab + (uufin(i)-estac(i+1))**2
enddo
eab= dsqrt(eab)
C guarda ultima solucion
if (ern.LE.distan) then
C write(13,30)n,cont,ele,kas,ern,eab,(cont*kas/60d0/60d0)
C write(*,30)n,cont,ele,kas,ern,eab,(cont*kas/60d0/60d0)
endif
C verifica lmin
if (ele.le.2.60d0)then
out1=.false.
endif
enddo
enddo
enddo
enddo
close(13)
20 format(6F12.6)
30 format(I4,I6,D14.7,D14.7,D14.7,D14.7,D14.7)
end
C
SUBROUTINE jacobi(a,n,np,d)
INTEGER n,np,nrot,NMAX
Double Precision a(np,np),d(np),v(np,np)
PARAMETER (NMAX=500)
INTEGER i,ip,iq,j
Double Precision c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
do 12 ip=1,n
do 11 iq=1,n
v(ip,iq)=0.0d0
11 continue
v(ip,ip)=1.0d0
12 continue
do 13 ip=1,n
b(ip)=a(ip,ip)
d(ip)=b(ip)
z(ip)=0.
13 continue
nrot=0
do 24 i=1,50
sm=0.
154
do 15 ip=1,n-1
do 14 iq=ip+1,n
sm=sm+dabs(a(ip,iq))
14 continue
15 continue
if(sm.eq.0.)return
if(i.lt.4)then
tresh=0.2d0*sm/n**2
else
tresh=0.0d0
endif
do 22 ip=1,n-1
do 21 iq=ip+1,n
g=100.*dabs(a(ip,iq))
if((i.gt.4).and.(dabs(d(ip))+
*g.eq.dabs(d(ip))).and.(dabs(d(iq))+g.eq.dabs(d(iq))))then
a(ip,iq)=0.
else if(dabs(a(ip,iq)).gt.tresh)then
h=d(iq)-d(ip)
if(dabs(h)+g.eq.dabs(h))then
t=a(ip,iq)/h
else
theta=0.5d0*h/a(ip,iq)
t=1.0d0/(dabs(theta)+dsqrt(1.0d0+theta**2))
if(theta.lt.0.)t=-t
endif
c=1.0d0/sqrt(1+t**2)
s=t*c
tau=s/(1.0d0+c)
h=t*a(ip,iq)
z(ip)=z(ip)-h
z(iq)=z(iq)+h
d(ip)=d(ip)-h
d(iq)=d(iq)+h
a(ip,iq)=0.
do 16 j=1,ip-1
g=a(j,ip)
h=a(j,iq)
a(j,ip)=g-s*(h+g*tau)
16 continue
do 17 j=ip+1,iq-1
155
g=a(ip,j)
h=a(j,iq)
17 continue
do 18 j=iq+1,n
g=a(ip,j)
h=a(iq,j)
a(iq,j)=h+s*(g-h*tau)
18 continue
do 19 j=1,n
g=v(j,ip)
h=v(j,iq)
v(j,ip)=g-s*(h+g*tau)
v(j,iq)=h+s*(g-h*tau)
19 continue
nrot=nrot+1
endif
21 continue
22 continue
do 23 ip=1,n
b(ip)=b(ip)+z(ip)
d(ip)=b(ip)
z(ip)=0.
23 continue
24 continue
pause ’too many iterations in jacobi’
return
END
C
SUBROUTINE eigsrt(d,n,np)
INTEGER n,np
Double Precision d(np),v(np,np)
INTEGER i,j,k
Double Precision p
do 13 i=1,n-1
k=i
p=d(i)
do 11 j=i+1,n
if(d(j).ge.p)then
156
k=j
p=d(j)
endif
11 continue
if(k.ne.i)then
d(k)=d(i)
d(i)=p
do 12 j=1,n
p=v(j,i)
v(j,i)=v(j,k)
v(j,k)=p
12 continue
endif
13 continue
return
END
C
SUBROUTINE choldc(a,n,np,p)
INTEGER n,np
Double Precision a(np,np),p(np)
INTEGER i,j,k
Double Precision sum
do 13 i=1,n
do 12 j=i,n
sum=a(i,j)
do 11 k=i-1,1,-1
sum=sum-a(i,k)*a(j,k)
11 continue
if(i.eq.j)then
p(i)=dsqrt(sum)
else
a(j,i)=sum/p(i)
endif
12 continue
13 continue
do 15 i=1,n
do 14 j=i+1,n
a(i,j)=a(j,i)
14 continue
a(i,i)=p(i)
15 continue
157
return
END
C
SUBROUTINE ernor(np,n,um1,um,er)
Double Precision um1(np),um(np),er,su1,su2
Integer i
su1= 0.0d0
su2= 0.0d0
do i=1,n,1
su1= su1 + (um1(i)-um(i))**2
enddo
er=dsqrt(su1)
END
C
SUBROUTINE subpaso(np,n,m1,m2,sant,vect,sact)
Integer n,np,i,j,k
Double Precision m1(np,np),m2(np,np),sant(np),vect(np),
* sact(np),b(n),c(n)
C calcula b = m2*sant + vect
do i=1,n,1
b(i)= vect(i)
do j=1,n,1
b(i)= b(i) + m2(i,j)*sant(j)
enddo
enddo
C resuelve c por triangular inferior: m1*c=b
do i=1,n,1
c(i)=b(i)
do k=1,i-1,1
c(i)=c(i)-m1(i,k)*c(k)
enddo
c(i)=c(i)/m1(i,i)
enddo
C resuelve sact por triangular superior: m1*sact=c
do i=1,np,1
sact(i)=0.0d0
enddo
do i=n,1,-1
sact(i)=c(i)
do k=i+1,n,1
sact(i)=sact(i)-m1(i,k)*sact(k)
enddo
158
sact(i)=sact(i)/m1(i,i)
enddo
return
END
C
TAREA 7. Modifique el programa anterior de tal forma que pueda utilizarse en una sola ejecu-
ción para discretizaciones n = 50, 100, 150, ..., 500. Entrege impreso y archivo electrónico (*.f) para
verificación.
2.11.2. Douglas − Rachf ord.

Bajo las mismas consideraciones dadas en el algoritmo semi-implícito de Euler, procedemos a la
simulación para el algoritmo de Douglas-Rachford.
0.Determine β y µ : β = A y µ = M .

 Seleccione ω tal que: 0 ≤ ω ≤ 1, o
f
 Determine ω como: ω = , donde u es solución de Au = f .
A u
Dado un error relativo normalizado1 límite ε suficientemente pequeño fijo.
 Seleccione ℓ tal que: 0 < ε < ℓ ≤ 2 (1 + ω) .


 Determine k tal que:
/
o

 −2 (ℓ − (1 + ω)) + 2 (1 + ω) ((1 + ω) − 2ℓ) µ

 k= .
 ℓ β
 Seleccione ℓ tal que: 0 < ε$< ℓ < (1 + ω) %.
ℓ µ (2.77)
 Determine k tal que: k = .
(1 + ω) − ℓ β
1. Dado u0 ,
2. Conocido um , para m ≥ 0,
3. Calcular
um+1 y um+1 ,
−1
um+1 = M + 12 kA M − 12 kA um + kf m+1 ,
−1
um+1 = M + 12 kA Mum − 12 kAum+1 + kf m+1 .
um+1 − um
4. Si ≤ ε, entonces terminar,
um

Además tenemos dos parametrizaciones:

 Seleccione ℓ tal que: 0 < ε$< ℓ < (1 + ω)
%.
ℓ µ
(1 + ω) − ℓ β
159
es llamada Douglas-Rachford, mientras que,


 Seleccione ℓ tal que: 0 < ε < ℓ ≤ 12 (1 + ω) .

/

 −2 (ℓ − (1 + ω)) + 2 (1 + ω) ((1 + ω) − 2ℓ) µ

 k= .
ℓ β
es llamada Douglas-Rachford2.
El comportamiento numérico de las triangulaciones del dominio n = 5, ..., 100, bajo las dos parame-
trizaciones es similar, hasta un error relativo normalizado del orden de ε = 1E − 12, en el caso de
la parametrización Douglas-Rachford2 tenemos que es insensible al valor de ω que se tome o selec-
cione. Además en ambas parametrizaciones tenemos que el valor mínimo de iteraciones para cada
discretización es igual.
Fig. 2.22
En la figura 2.22 tenemos la simulación de la triangulación n = 100, ω = 0, ε = 1E −8, con parámetro

l = 0.045 equivalente a un tamaño de paso en el tiempo k = 159.8550 seg.
Comparando las triangulaciones n=100, n=50 y n=5, tenemos la siguiente tabla de valores:
160
w=0 n = 100 n = 50 n=5
ε iter iter iter
1E − 2 230 227.8908 14.5597 132 411.2660 15.0798 22 3179.880 19.4326
1E − 3 351 198.5579 19.3594 193 365.6536 19.6031 29 3101.0740 24.9809
1E − 4 474 182.8369 24.0735 255 344.7270 24.4182 36 3072.3240 30.7232
1E − 5 599 172.5746 28.7145 317 335.1197 29.5092 43 3053.3930 36.4711
1E − 6 724 166.3440 33.4536 380 326.0226 34.4135 49 3169.8540 43.1452
1E − 7 849 163.3926 38.5334 444 317.3960 39.1455 56 3130.2550 48.6929
1E − 8 975 160.5430 43.4804 507 314.6188 44.3088 63 3110.7520 54.4382
1E − 9 1101 157.7900 48.2574 570 311.8887 49.3824 70 3091.4430 60.1114
1E − 10 1227 155.1288 52.8731 634 306.5649 53.9895 76 3169.8540 66.9191
1E − 11 1354 155.1288 58.3457 698 306.5649 59.4395 83 3130.2550 72.1698
1E − 12 1502 152.5548 63.6492 764 303.9690 64.5090 90 3120.4800 78.0120
Tanto para la parametrización Douglas-Rachford como Douglas-Rachford2, se tiene que el tiempo de

llegada al estado estacionario para n = 5 está retrasado.
2.11.3. P eaceman − Rachf ord.

Bajo las mismas consideraciones dadas en los algoritmos anteriores, procedemos a la simulación para
el algoritmo de Peaceman-Rachford.

f
A u
Dado un error relativo normalizado límite ε suficientemente pequeño fijo.

 Seleccione ℓ tal que: 0 < ε < ℓ ≤ (1 + ω) .

/
o

 −4 (ℓ − 2 (1 + ω)) + 8 (1 + ω) ((1 + ω) − ℓ) µ

 k= .
 ℓ β (2.78)
 Seleccione ℓ tal que: 0 < ε$< ℓ < 2 (1 + ω) %.
2ℓ µ
2 (1 + ω) − ℓ β
1. Dado u0 ,
1
3. Calcular
 um+ 2 y um+1 ,

 um+ 12 = M + 1 kA −1 M − 1 kA um + 1 kf m+ 12 ,
4 4 2

 um+1 = M + 1 kA −1 M − 1 kA um+ 12 + 1 kf m+1 .
4 4 2
161
um+1 − um
um


 Seleccione ℓ tal que: 0 < ε$< ℓ < 2 (1 + ω)
%.
2ℓ µ
2 (1 + ω) − ℓ β
es llamada Peaceman-Rachford, mientras que,


 Seleccione ℓ tal que: 0 < ε < ℓ ≤ (1 + ω) .

/

 −4 (ℓ − 2 (1 + ω)) + 8 (1 + ω) ((1 + ω) − ℓ) µ

 k= .
ℓ β
es llamada Peaceman-Rachford2.
El comportamiento numérico de las triangulaciones del dominio n = 5, ..., 100, bajo las dos para-
metrizaciones es similar, hasta un error relativo normalizado del orden de ε = 1E − 12, en el caso
de la parametrización Peaceman-Rachford2 tenemos que es insensible al valor de ω que se tome o
seleccione. Además en ambas parametrizaciones tenemos que el valor mínimo de iteraciones para
cada discretización es igual.
Fig. 2.23
162
En la figura 2.23 tenemos la simulación de la triangulación n=100, ω = 0, ε = 1E − 8, con parámetro
w=0 n = 100 n = 50 n=5

ε iter iter iter
1E − 2 132 428.1247 15.6979 74 781.9958 16.0744 11 5624.7250 17.1867
1E − 3 193 376.2885 20.1733 104 721.6510 20.8477 13 6175.4260 22.3002
1E − 4 254 352.4056 24.8642 135 678.3253 25.4372 16 6006.7690 26.6968
1E − 5 316 339.4239 29.7939 166 651.9999 30.0644 19 5964.5410 31.4795
1E − 6 378 327.3269 34.3693 197 639.5244 34.9962 21 6217.5550 36.2691
1E − 7 440 323.4757 39.5359 228 627.4756 39.7401 24 6133.2880 40.8886
1E − 8 502 316.0270 44.0682 259 619.6694 44.5818 27 6048.9650 45.3672
1E − 9 565 312.4241 49.0332 290 615.8319 49.6087 29 6259.6820 50.4252
1E − 10 627 308.8986 53.7998 321 612.0372 54.5733 32 6175.4260 54.8927
1E − 11 690 308.8986 59.2056 352 604.5732 59.1138 35 6133.2880 59.6292
1E − 12 823 274.5899 62.7743 384 600.9026 64.0963 37 6259.6820 64.3356
Tanto para la parametrización Peaceman-Rachford como Peaceman-Rachford2, se tiene que el tiempo

de llegada al estado estacionario para las triangulaciones y ordenes de error relativo normalizado son
similares.
2.11.4. Esquema − theta.

Bajo las mismas consideraciones dadas en los algoritmos anteriores, procedemos a la simulación para
el algoritmo θ-Glowinski.

f
A u
Dado un error relativo normalizado límite ε suficientemente pequeño fijo.
 √

 Seleccione ℓ tal que: 0 < ε < 3 (1 + ω) < ℓ <
8 4 79−26
3
(1 + ω) .


 Determine
$ k tal que: %
0
2 2 o (2.79)

 − (3ℓ − 14 (1 + ω)) + −3ℓ − 52ℓ (1 + ω) + 196 (1 + ω)

 µ
 k=3
 (3ℓ − 8 (1 + ω)) β
 Seleccione ℓ tal que: 0 < ε$< ℓ < 3 (1 + ω)%
14
.
6ℓ µ
14 (1 + ω) − 3ℓ β
1. Dado u0 ,
163
1 2
3.Calcular um+ 3 , um+ 3 y um+1 ,
$ %−1 $$ % %

 1 k k k m+ 1

 m+
u 3 = M+ A m
M− A u + f 3 ,

 6 6 % 3

 $ %−1 $$ %
m+ 23 k k m+ 13 k m+ 2
u = M+ A M− A u + f 3 ,

 6 6 3

 $ %−1 $$ % %

 k k m+ 23 k m+1

 u m+1
= M+ A M− A u + f .
6 6 3
um+1 − um
um


 Seleccione ℓ tal que: 0 < ε$< ℓ < 14
3
(1 + ω)%
.
6ℓ µ
14 (1 + ω) − 3ℓ β
es llamada Theta, mientras que,
 √

 Seleccione ℓ tal que: 0 < ε < 8
(1 + ω) < ℓ ≤ 4 79−26
(1 + ω) .

 3 3

$ 0 %
 2 2

 − (3ℓ − 14 (1 + ω)) + −3ℓ − 52ℓ (1 + ω) + 196 (1 + ω)

 µ
 k=3
(3ℓ − 8 (1 + ω)) β
es llamada Theta2.
El comportamiento numérico de las triangulaciones del dominio n = 5, ..., 100, bajo las dos para-
metrizaciones es similar, hasta un error relativo normalizado del orden de ε = 1E − 12. Además
en ambas parametrizaciones tenemos que el valor mínimo de iteraciones para cada discretización es
igual.
164
Fig. 2.24
En la figura 2.24 tenemos la simulación de la triangulación n=100, ω = 0, ε = 1E − 8, con parámetro
w=0 n = 100 n = 50 n=5

ε iter iter iter
1E − 2 131 418.282 15.221 67 879.319 16.365 8 8323.632 18.497
1E − 3 172 418.282 19.985 88 857.144 20.952 10 8323.632 23.121
1E − 4 213 418.282 24.748 108 857.144 25.714 11 9394.302 28.705
1E − 5 255 418.282 29.628 129 857.144 30.714 13 9074.306 32.768
1E − 6 297 418.282 34.508 149 857.144 35.476 15 8951.482 37.298
1E − 7 339 418.282 39.388 170 857.144 40.476 17 8891.127 41.986
1E − 8 380 418.282 44.152 191 857.144 45.476 19 8831.462 46.610
1E − 9 422 418.282 49.032 211 857.144 50.238 20 9328.772 51.827
1E − 10 464 418.282 53.912 232 857.144 55.238 22 9200.035 56.222
1E − 11 507 418.282 58.908 254 857.144 60.476 24 9136.801 60.912
1E − 12 1587 140.770 62.056 289 836.042 67.116 26 9074.306 65.537
Tanto para la parametrización Theta como Theta2, se tiene que el tiempo de llegada al estado
estacionario para las triangulaciones y ordenes de error relativo normalizado son similares.
165

MEF035A165

Cargado por

Información del documentohacer clic para expandir la información del documento

Copyright:

Formatos disponibles

MEF035A165

Cargado por

Información del documento

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

MEF035A165

Cargado por

Copyright:

Formatos disponibles

2.

Modelación y Simulación Numérica de Problemas de Conducción de

Figura 1. Modelación de Problemas.

Figura 2. Esquema Conceptual.

Fig. 2.1. Problema de Valores sobre la Frontera.

Gráfica 1. Balance de Calor en el Interior.

Gráfica 2. Temperatura Prescrita.

Gráfica 3. Flujo Prescrito.

2.2. Modelo físico

Identificamos los términos del modelo como:

u (x, 0) = u0 (x) , x ∈ Ω, (2.5)

Notar en gráficas (4 − 6) y en el modelo (2.5 − 2.8) se ha cambiado la temperatura u por la temprea-

2.3. Funciones primitivas.

Gráfica de la función p. Gráfica de la función Dp.

Gráfica de la función Dp. Gráfica de la función Dp dx.

En el interior del cuerpo.

Gráfica 4a. Gráfica 4b.

En la frontera de temperatura prescrita.

Grafica 5a. Gráfica 5b.

En la frontera de flujo prescrito.

Y los correspondientes subpotenciales φ, ψ 1 y ψ 2 están dados por las expresiones siguientes:

φ (x, t; u (x, t)) = 0, u (x, t) > 0, (x, t) ∈ Ω × [0, T ] (2.9)

2.4. Modelos Matemáticos.

∂f (u) = u∗ ∈ E ′ : f (w) ≥ f (u) + u∗ , w − uE ′ ×E ∀w ∈ E . (2.12)

Observación. f : E −→ R ∪ {+∞} es convexa si y sólo si su subdiferencial ∂f : R −→ 2R es

u∗ ∈ ∂f (u) ⊂ E ′ ⇐⇒ f (w) ≥ f (u) + u∗ , w − uE′ ×E ∀w ∈ E

De lo expresado en la observación anterior, procedemos a modificar el modelo (2.5 − 2.8) y expresamos

γu (x, t) ∈ D (ψ 1 (x, t; ·)) = {û (x, t)} :

γu (x, t) ∈ D (ψ 1 (x, t; ·)) = {û (x, t)} :

Encuentre u tal que se cumple:

u (t)∈ D (Φ (t; ·)) :

 ∈ D (Ψ1 (t; ·)) = {û (t)} :

γu (t) ∈ D (Ψ2 (t; ·)) :

{(−K grad u (t)) · n} {γη (t) − γu (t)} dΓ, ∀η (t) ∈ Kt .

aplicando teorema de la divergencia al término izquierdo tenemos:

sean q = −K grad u y α = {η − u} obtenemos la fórmula de Green correspondiente

Aplicando la fórmula de Green [2.31] y simplificando obtenemos la formulación variacional global

2.5. Modelo Discreto.

Dada la condición inicial discreta:

Substituyendo y simplificando obtenemos la versión coordenada del modelo discreto.

Dada la condición inicial discreta:

Ktmh = η : R → Rmh | η (t) > 0 en Ωh , γη (t)|Γ1h = û sobre Γ1h . (2.40)

Dada la condición inicial discreta:

Ktmh = η : R → Rmh | η (t) > 0 en Ωh , γη (t)|Γ1h = û sobre Γ1h . (2.44)

Ktmh = η : R → Rmh | η (t) > 0 en Ωh , γη (t)|Γ1h = û sobre Γ1h . (2.47)

Dada la condición inicial discreta:

Ktmh = η : R → Rmh | η (t) > 0 en Ωh , γη (t)|Γ1h = û sobre Γ1h . (2.50)

De acuerdo a la versión coordenada modelo semi-discreto espacial:

y donde expresamos la combinación convexa de A, con α ∈ (0, 1) ,

El método de Galerkin, para la aproximación del problema (2.53)

En particular para éste propósito estudiaremos técnicas de elemento finito.

Una vez establecida la triangulación τ h sobre Ω, se define

a este espacio se le conoce como Espacio de Elemento Finito Global .

estos espacios son llamados Espacios

-Por Elemento Finito entendemos a la terna

Caso lineal escalar.

el espacio de elemento finito local PE es

y el conjunto de grados de libertad ΣE es

ΣE = {φi ∈ L (PE , R) : φ1 (p) ∈ PE , φ2 (p) ∈ PE }

∂f (u) = u∗ ∈ E ′ : f (w) ≥ f (u) + u∗ , w − uE ′ ×E ∀w ∈ E . (2.12)

u∗ ∈ ∂f (u) ⊂ E ′ ⇐⇒ f (w) ≥ f (u) + u∗ , w − uE′ ×E ∀w ∈ E

u (t)∈ D (Φ (t; ·)) :

∈ D (Ψ1 (t; ·)) = {û (t)} :

El método de Galerkin, para la aproximación del problema (2.53)

para la primera función base tenemos:

con los valores de las propiedades físicas y dimensiones siguientes: