Analisis-y-solucion-de-la-ecuacion-de-momentum-y-energa-utilizando-tecnicas-anlticas-y-el-metodo-de-elemento-finito

Vista previa del material en texto

UNIVERSIDAD NACIONAL AUTÓNOMA DE 
MÉXICO 
FACULTAD DE QUÍMICA 
 
 ANÁLISIS Y SOLUCIÓN DE LA ECUACIÓN DE DIFUSIÓN DE MOMENTUM Y ENERGÍA 
USANDO TÉCNICAS ANALÍTICAS Y EL MÉTODO DE ELEMENTO FINITO. 
 
T E S I S 
 
QUE PARA OBTENER EL TÍTULO DE: 
INGENIERO QUÍMICO 
 
PRESENTA: 
FERNÁNDEZ MARBÁN RAÚL 
 
CIUDAD DE MÉXICO 2016
 
UNAM – Dirección General de Bibliotecas 
Tesis Digitales 
Restricciones de uso 
 
DERECHOS RESERVADOS © 
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL 
 
Todo el material contenido en esta tesis esta protegido por la Ley Federal 
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). 
El uso de imágenes, fragmentos de videos, y demás material que sea 
objeto de protección de los derechos de autor, será exclusivamente para 
fines educativos e informativos y deberá citar la fuente donde la obtuvo 
mencionando el autor o autores. Cualquier uso distinto como el lucro, 
reproducción, edición o modificación, será perseguido y sancionado por el 
respectivo titular de los Derechos de Autor. 
 
 
 
2 
 
JURADO ASIGNADO: 
PRESIDENTE: MANUEL VAZQUEZ ISLAS 
VOCAL: ALBERTO CASTELLANOS CAMPILLO 
SECRETARIO: LEONARDO DAMIAN SORIA RODRIGUEZ 
1er. SUPLENTE: ILEANA RODRIGUEZ CASTAÑEDA 
2° SUPLENTE: AMPARO MAYA ROMERO 
 
 
SITIO DONDE SE DESARROLLÓ EL TEMA: DEPARTAMENTO DE INGENIERIA QUIMICA - 
CIUDAD UNIVERSITARIA – CIUDAD DE MEXICO 
 
 
ASESOR DEL TEMA: LEONARDO DAMIAN SORIA RODRIGUEZ 
 
Firma 
 
 
SUSTENTANTE: RAÚL FERNÁNDEZ MARBÁN 
 
Firma 
 
 
 
 
 
 
3 
 
Índice general 
 
1. Introducción ........................................................................................................ 8 
1.1. Ecuaciones diferenciales parciales ................................................................ 8 
1.2. El orden de una EPD...................................................................................... 8 
1.3. Teoremas para la solución de ecuaciones diferenciales .............................. 10 
1.4. Condiciones a la frontera ............................................................................. 10 
1.4.1 Condiciones a la frontera de primera clase o Dirichlet ........................... 10 
1.4.2 Condiciones a la frontera de segunda clase o Neumann ...................... 11 
1.5 Método de separación de variables.............................................................. 11 
2. Antecedentes .................................................................................................... 13 
2.1 Funciones ortogonales ................................................................................. 13 
2.2 Series de Fourier .......................................................................................... 14 
2.3 Teorema de Sturm Liouville ......................................................................... 15 
2.4 Ecuaciones de Navier-Stokes ...................................................................... 17 
2.5 Antecedentes - Método de elemento finito (FEM) ........................................ 18 
2.6 Independencia lineal .................................................................................... 18 
2.7 Métodos variacionales de Rayleight Ritz ...................................................... 19 
2.8 Mallado ......................................................................................................... 20 
3 Objetivos ........................................................................................................... 22 
3.1 Objetivos generales...................................................................................... 22 
3.2 Objetivos particulares ................................................................................... 22 
4 Desarrollo del método analítico ...................................................................... 23 
4.1 Caso I Ecuación de difusión ......................................................................... 23 
4.2 Caso II – Ecuación de transferencia de momentum, con cambio en las C.F25 
4.3 Caso III – Ecuación de difusión de energía para el caso de 2 dimensiones 27 
4.4 Caso IV – Ecuación de difusión de energía para el caso de 3 dimensiones 29 
4.5 Caso V – Ecuación de difusión de masa ...................................................... 32 
4.6 Caso VI – Método perturbativo y capa límite ................................................ 35 
5 Método del elemento finito (MEF ó FEM) ....................................................... 38 
5.1 Solución débil ............................................................................................... 38 
5.2 Acercamiento de Galerkin ............................................................................ 41 
5.3 Matriz de elementos (Element stiffness matrix) para el caso 1D .................. 44 
4 
 
5.4 Uso Matlab para el caso numérico el caso 1D (Elemento barra) ................. 47 
5.5 Caso 2 D (Mesh triangular) .......................................................................... 51 
5.6 Ecuación diferencial parcial 2D .................................................................... 56 
6 Creación del mallado (Mesh) ........................................................................... 58 
6.1 Bresenham ................................................................................................... 58 
6.2 Set Domian .................................................................................................. 60 
6.3 Flood fill ........................................................................................................ 61 
6.4 Create triangles ............................................................................................ 63 
6.5 Create Mesh ................................................................................................. 64 
6.6 Caso 2D (Mesh cuadrado) ........................................................................... 66 
6.7 Cuadratura Gaussiana ................................................................................. 68 
7 Comparación del método analítico caso IV – Ecuación de transferencia de 
energía - Contra el método numérico de elemento (FEM) con diferencias finitas 
de Euler .................................................................................................................... 69 
7.1 Diferencias finitas ......................................................................................... 70 
8 Conclusiones y recomendaciones .................................................................. 75 
9 Trabajo futuro ................................................................................................... 77 
10 Apéndice ............................................................................................................ 78 
10.1 Implementación de los sistemas de cómputo, herramienta Matlab para el 
método de elemento finito. ..................................................................................... 78 
10.2 Caso 1D (Mallado de barras) ....................................................................... 78 
10.2 Mallado caso 2D ............................................................................................ 80 
10.3 Caso 2D mallado cuadrático ........................................................................ 87 
11 Bibliografía y referencias APA .......................................................................... 91 
 
 
 
 
 
5 
 
Tabla de figuras 
 
Figura 1 - Series de Fourier .................................................................................... 14 
Figura 2 - Nodos ...................................................................................................... 21 
Figura 3 - Mallado .................................................................................................... 21 
Figura 4 – Transición ..............................................................................................21 
Figura 5 - Placas ...................................................................................................... 32 
Figura 6 - Banda elástica ........................................................................................ 39 
Figura 7 - Metodología FEM .................................................................................... 41 
Figura 8- Barra Segmentada ................................................................................... 42 
Figura 9- División de secciones ............................................................................. 42 
Figura 10- Funciones de forma .............................................................................. 43 
Figura 11- Solución desplazamiento ..................................................................... 43 
Figura 12- Matriz FEM ............................................................................................. 45 
Figura 13- No. grados de libertad .......................................................................... 45 
Figura 14- Comparación sol. analítica vs numérica ............................................. 51 
Figura 15 -Elemento dominio ................................................................................. 52 
Figura 16 – Elementos triangulares ....................................................................... 54 
Figure 17 - Mapping table figura elemento izquierdo ........................................... 54 
Figure 18 - Mapping Table elemento derecho ....................................................... 54 
Figura 19 - Matriz global de elementos.................................................................. 55 
Figura 20 - Solución perfil de Temperaturas FEM ................................................ 58 
Figura 21 - Bresenham ............................................................................................ 58 
Figure 22 - Desplazamineto de puntos .................................................................. 60 
Figura 23 - Set Domain - Bresenham ..................................................................... 61 
Figura 24 – Floodfill y Seed point .......................................................................... 62 
Figura 25 - Orientación de triángulos .................................................................... 63 
Figura 26 - Selección de puntos ............................................................................ 64 
Figura 27 - Trazado de elemento ............................................................................ 65 
Figura 28 - Llenado Floodfill ................................................................................... 65 
Figura 29 - Elemento con límites Dirichlet ............................................................ 65 
Figura 30 - Cambio de coordenadas ...................................................................... 66 
Figure 31 - Barra expuesta a una fuente externa de calor ................................... 70 
Figura 32 - Primer plano - Diferencias finitas ....................................................... 74 
Figura 33 - Segundo plano - Diferencias finitas ................................................... 74 
Figura 34 - Tercer plano - Diferencias finitas ........................................................ 74 
Figura 35 - Equilibrio térmico al tiempo 33 ........................................................... 75 
 
 
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630312
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630322
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630323
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630331
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630332
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630333
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630334
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630337
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630338
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630339
file:///C:/Users/FERNARAK/Desktop/Tesis%20III.docx%23_Toc452630340
6 
 
Resumen 
Éste trabajo, titulado “Análisis y solución de la ecuación de disfusión de momentum y 
energía usando técnicas analíticas y el método de elemento finito’’ propone una serie de 
estrategias para facilitar el cálculo analítico y aproximación numérica a las ecuaciones 
diferenciales parciales. Se presentará en dos partes, una la parte analítica que produce 
soluciones exactas pero se requiere de más abstracción y cálculo y otra que será una 
aproximación numérica (Método del elemento finito - MEF) a la solución con apoyo del 
software MatLab. Ambas métodos si bien, son muy usados presentan sus ventajas uno 
con respecto al otro. La computadora tomará un rol importante en éste trabajo, pues se 
requiere de cálculos iterativos respecto de ecuaciones diferenciales parciales de segundo 
orden o superior. 
El trabajo será de alto impacto, pues las ecuaciones que se analizarán, han demostrado 
tener muy buenos resultados en el campo de la reología, si bien es cierto que éste 
proyecto apunta más a la matemática teórica, porque tomará conceptos analizados y 
estudiados a lo largo de la carrera, como lo son: el de la capa limite, perfil de velocidades, 
y los procesos difusivos, también lo es que en la parte analítica, se necesita considerar el 
cálculo por aproximaciones numéricas o diferencias finitas, logrado también a la par con el 
método del elemento finito. La técnica (MEF) se propone para conseguir el objetivo final 
de ésta tesis y será su punto medular, pues pretende ser una perspectiva diferente para 
abordar la solución de ecuaciones diferenciales, que podría ser de utilidad para su 
implementación en MatLab; es de subrayarse, que dicha técnica, nos remonta a 
conceptos básicos de algebra lineal, entre ellos espacios vectoriales, euclidianos, 
combinación lineal, hasta la discretización de las ecuaciones diferenciales parciales en 
elementos más sencillos, es decir que puedan ser descompuestas en sistemas de 
ecuaciones algebraicas, para posteriormente resolver el sistema con el uso de matrices. 
En el mundo actual se cuenta con muchas herramientas, softwares, calculadoras y 
métodos analíticos que dan solución a la mayoría de los problemas cotidianos, sin 
embargo, siempre habrá excepciones y para ellas existen los métodos numéricos, que 
gracias a su gran poder y el uso de condiciones a la frontera (C.F) se pueden aproximar 
las soluciones con un error asociado, pero despreciable a ecuaciones de mayor grado de 
dificultad o aquellas para las cuales no se conocen métodos de solución de forma 
desarrollada o analítica. 
Todos estos procedimientos y conocimientos previos y recientes de matemáticas serán 
aplicados a los modelos o ecuaciones diferenciales de momentum, las cuales describen 
fluidos viscosos o de viscosidad nula, difusión, esfuerzos cortantes, entre muchos otros 
factores que ayudan a comprender la mecánica de fluidos, que es determinante en la 
industria química. 
 
 
7 
 
Agradecimientos 
Primero que a nadie me gustaría agradecer sinceramente a mi asesor de tesis el Ing. 
Químico Leonardo Damián Soria Rodríguez, por su apoyo incondicional e interés para el 
desarrollo de éste trabajo. Su orientación y conocimiento, me permitieron tomar rumbo y 
decisión para encaminar mis pasos en la investigación que sustenta éste material. a 
través del cual deseo titularme. Desde que aceptó la solicitud para tomar su curso a 
distancia, durante mi estancia académica en el extranjero, ha demostrado una disposición 
absoluta como guía y mentor, razón por la que mi gratitud y admiración es permanente. 
Gracias también al profesor Manuel Vázquez Islas, a quien tengo en muy alta estima, por 
su apoyo y amistad a lo largo de mi formación, pues a pesar de no haber tomado ninguna 
materia con él, siempre se presentó abierto a ayudarme en asignaturas como cálculo IIy 
ecuaciones diferenciales. Le reconozco asimismo la orientación en la elaboración del 
manual “Apuntes MOI 1,2,3 resolución de ecuaciones diferenciales” como medio de 
liberación del servicio social. 
Agradezco además al Ing. Alberto Castellanos Campillo y al jurado designado, por 
tomarse el tiempo para la revisión de éste trabajo y formar parte del jurado evaluador del 
examen profesional que corresponde en el caso. 
De manera especial dedico éste proyecto a mi madre María Eugenia Marbán, porque fue 
y es el principal pilar en mi vida y formación profesional, sentó en mí las bases de 
responsabilidad, honestidad y cariño, al heredarme sus deseos de superación personal, 
así como también sus virtudes infinitas que me harán estar por siempre agradecido. 
A mi padre Eduardo H. Fernández y a mis hermanos Eduardo y Cristina Fernández 
Marbán, por haber formado en parte mi carácter y por estar siempre ahí para mí. 
Y por último y gracias a ésta máxima casa de estudios UNAM, responsable en buena 
medida de mi formación académica, incluyendo el intercambio académico internacional 
tan formativo, en el ámbito profesional. Estoy orgulloso de haber formado parte de su 
alumnado. 
 
 
 
 
 
8 
 
1. Introducción 
1.1. Ecuaciones diferenciales parciales 
 
La comprensión de los fenómenos de la naturaleza, requiere de la ayuda de modelos 
matemáticos, es decir de ecuaciones diferenciales tanto ordinarias como parciales, para 
las áreas de física, matemáticas e ingeniería. Así pues en éste trabajo estudiaremos 
diversas ecuaciones diferenciales parciales, pero en particular, aquellas que tienen una 
aplicación directa en la Ing. Química. 
Dentro de éstas ecuaciones diferenciales parciales (EDP) tenemos las elípticas, 
hiperbólicas y parabólicas, para resolverlas, se requerirá del uso de algunas técnicas 
importantes, como son: Separación de variables, simetrizaciones, conjuntos lineales, 
solución de matrices y los métodos variacionales de Galerkin a su vez para el método 
numérico de elemento finito (MEF) para éste último, se propone una técnica no muy 
socorrida en la facultad de química, que consiste en la discretización de las ecuaciones 
diferenciales parciales a algebraicas. 
Para que el lector se familiarice, con la notación usada en el siguiente trabajo, respecto a 
las ecuaciones diferenciales parciales, es pertinente, considerar los siguientes conceptos: 
¿Qué es una ecuación diferencial parcial? 
Es un modelo matemático que contiene derivadas parciales, cuyas incógnitas a estudiar, 
dependen de diversas variables, es decir en la ecuación estarán presentes las funciones 
propias y sus derivadas. Suelen ser ecuaciones que involucran funciones U que 
dependen de variables independientes: x, y, z, t, así como las derivadas de U. Se usan 
con mayor frecuencia en: problemas de electromagnetismo, transferencia de calor, 
propagación de calor, membranas plásticas, electromecánica, magnetismo, campos 
gravitacionales, entre otras muchas aplicaciones. Estas ecuaciones presentan ésta forma: 
F(x1, x2, x3, … , xn, n,
∂u
∂x1
, … ,
∂u
∂xn
,
∂un
n
∂
k1
x1 ∂
k2
x2 …∂
kn
xn
) = 0 
1.2. El orden de una EDP 
 
Se le considera a aquel exponente de la derivada más grande, que figura en la ecuación. 
Como ejemplo se puede proponer una de segundo orden, con forma: 
∂
2
U
∂x
2 − 
∂
2
U
∂y
2 = 0 
 
9 
 
Es importante resaltar, que una ecuación diferencial siempre da una familia de soluciones 
a diferencia de las algebraicas donde se obtienen resultados numéricos. 
Dentro de éstas EDP hay dos casos, una ecuación diferencial parcial, homogénea: 
(
∂
2
U
∂x
2)
2
+ (
∂
2
U
∂y
2)
2
= 0 ….. Presenta familias de soluciones reales 
Y otra que por el contrario no es homogénea sino que: 
(
∂
2
U
∂x
2)
2
+ (
∂
2
U
∂y
2)
2
+ 1 = 0 ….. Presenta familias de soluciones imaginarias. 
Para órdenes superiores al segundo, no se tienen soluciones analíticas, sino más bien 
numéricas, éste punto será abordado más adelante. Sin embargo no es sólo importante 
determinar el orden de la EDP, sino también el tipo de ecuación con la que se trabaja, es 
decir si se trata de ecuaciones: Hiperbólicas, parabólicas o elípticas (Sólo para 
ecuaciones diferenciales parciales de segundo orden, mayor se complica la geometría del 
sistema por tener polígonos irregulares, así como también deben de contar con 
coeficientes constantes y no variables) 
Para determinar el orden del tipo de ecuación que estamos tratando, será fundamental el 
espacio o dominio omega (Ω) o sea un dominio abierto en Rn que delimite las fronteras del 
sistema. 
Para un dominio Omega (plano x,y) con ciertas condiciones en el contorno o a la frontera, 
se dice que la ecuación es: 
Considerando (Λ) como el determinante del sistema. 
Hiperbólica en Ω si (Λ=B2-AC>0 en Ω) 
Ejemplo: 
∂2u
∂x2
− 9 
∂2u
∂x∂y
= 0 Donde A=1, B= -9 y C=0, se tiene que el Λ = 81>0, por ende se trata 
de una ecuación diferencial hiperbólica. 
Parabólica en Ω si (Λ=B2-AC< 0 en Ω) 
Ejemplo: 
∂2u
∂x2
+ 6 
∂2u
∂x∂y
+ 9 
∂2u
∂y2
= 0 Donde A=1, B= 6 y C=9, se tiene que el Λ = 0= 0, por 
ende se trata de una ecuación diferencial parabólica. 
 
 
10 
 
Elíptica en Ω si (Λ=B2-AC= 0 en Ω) 
Ejemplo: 
∂2u
∂x2
+
∂2u
∂x∂y
+ 
∂2u
∂y2
= 0 Donde A=1, B= 1 y C=1, se tiene que el Λ = -3< 0, por tanto se 
trata de una ecuación diferencial elíptica y no cualquiera, sino Laplace, abordada para el 
problema de elemento finito (MEF) 
1.3. Teoremas para la solución de ecuaciones diferenciales 
 
Considere L[u] = 0 , donde L es el operador diferencial. 
 Primer teorema: Si u(x,y) es la solución de la EDP L[u] = 0, entonces K*u, o sea un 
escalar(K) es también solución de la homogénea. 
 Segundo teorema: Si u(x,y) y u2(x,y) son solución de la EDP L[u] = 0, entonces 
u+u2, es también solución de la homogénea. 
 Tercer teorema (Principio de superposición): Es una herramienta matemática que 
permite descomponer un problema lineal en dos o más subproblemas más 
sencillos. 
1.4. Condiciones a la frontera 
Las ecuaciones diferenciales si bien representan muchos fenómenos fisicoquímicos, 
requieren de ciertas limitantes para soluciones específicas, éstas son denominadas 
condiciones a la frontera (también llamados como problemas de valor o de borde) Una 
solución de un problema de condiciones de frontera es una solución de una ecuación 
diferencial que también satisface condiciones de frontera. 
Problemas por ejemplo que involucran la ecuación de onda, son comúnmente problemas 
de condiciones de frontera, algunos otros importantes son los de Sturm-Liouville, tema 
que se explicará más adelante, donde su análisis necesita funciones propias y operadores 
diferenciales. 
Muchos de los primeros problemas de valor de frontera han sido estudiados mediante 
los problemas de Dirichlet, o buscando una función armónica (Solución de una ecuación 
de Laplace) cuya solución está dada por el principio de Dirichlet, así como también las 
condiciones a la frontera para las derivadas de Neumann. 
1.4.1 Condiciones a la frontera de primera clase o Dirichlet 
Se trata de una condición de frontera o contorno, denominado así en honor a Johann 
Peter Gustav Lejeune Dirichlet, cuando en una ecuación diferencial ordinaria o una 
en derivadas parciales, se le especifican los valores de la solución que necesita la frontera 
del dominio. 
https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_en_derivadas_parciales
11 
 
Casos especiales son: 
𝑇 = 𝑓(𝑟𝑠), 𝑇 = 𝑓(𝑡) ó 𝑇 = 𝑐𝑡𝑒 
1.4.2 Condiciones a la frontera de segunda clase o Neumann 
Son aquellas condiciones de frontera, también conocidas como condiciones contorno, 
llamadas así en alusión a Carl Neumann. Se presentan cuando a una ecuación 
diferencial ordinaria o en derivadas parciales, se le especifican los valores de 
la derivada de una solución tomada sobre la frontera o contorno del dominio. 
 
Casos especiales son: 
∂T
∂n
= f(rs) ; 
∂T
∂n
= f(t) ó 
∂T
∂n
=cte 
1.5 Método de separación de variables 
Es un procedimiento, para encontrar una solución particular para algunos problemas que 
involucran ecuaciones en derivadas parciales, como serie cuyos términos son el producto 
de funciones que presentan variables que pueden ser separadas por sencillos despejes. 
Uno de los métodos usados en la física y matemática para buscar soluciones a problemas 
descritos mediante ecuaciones diferenciales de derivadas parciales. Para las ecuaciones 
diferenciales parciales, funciona de manera diferente, es decir, para ecuaciones 
ordinarias, se pueden despejar funciones como producto o división, pero para las EDP, se 
busca una solución que sea un producto de funciones dependientes cada una de una sola 
variable. 
Para los casos hiperbólico y parabólico se buscará una solución de la forma: 
𝑈(𝑥, 𝑡) = 𝑀(𝑥)𝑁(𝑡) 
Procederemos a realizar un ejemplo con una ecuación muy concida en el campo de la 
ingenieria química: 
La ecuación de difusión presenta la siguiente forma 
 
∂T
∂t
= α2
∂2T
∂t2
 
 
Suponiendo las condiciones a la frontera: 
𝑇(0, 𝑡) = 0 ; 𝑇(1, 𝑡) = 0 ; 𝑇(𝑥, 0) = 𝛷(𝑥) 
 
12 
 
Como primer paso, se propone una función: 𝑇(𝑥, 𝑡) = 𝑀(𝑥)𝑁(𝑡) 
∂T
∂t
= M(x)N′(t) 
∂2T
∂t2
= M(x)N′′(t) 
Sustituyendo en la ecuación diferencial original de difusión: 
M(x)N′(t) = α2M(x)N′′(t) 
M′(x)
M(x)
= α2
N′′(t)
N(t)
= −𝜆2 
Aquí se aplicó el teorema de Sturm Loiuvielle. 
Se originan dos ecuaciones, hay que reconocer, que el número de ecuaciones 
generadas, será proporcional al grado de la EDP, dado que se trata de una 
ecuación de segundo orden, se tendrán dos ecuaciones: 
 (*) M’(x) + λ2M(x) = 0 y ($) α2N′′(t) + λ2N(t) = 0 
Como solución de la primera ecuación (*) se tiene 𝐌(𝐱) = 𝐀𝐞−𝛌𝟐𝐭 
Para la segunda ($) sólo se tiene que desarrollar un poco de algebra, y se resuelve 
como si se tratase de una ecuación cuadrática, es decir: 
N′′(t) +
λ2
α2
N(t) = 0 generá una nueva constante N′′(t) + kN(t) = 0 
Solución ($) N(t) = Bsen(xk) + Ccos(xk) valuando con C.F. 
Solución general, será el producto de las funciones: 𝐓(𝐱, 𝐭) = 𝐀𝐧𝐞−𝛌
𝟐𝐭𝐬𝐞𝐧(𝐧𝛑𝐱) 
Los temas expuestos anteriormente, presentan parte de la introducción y bases para la 
tesis a desarrollar, sin embargo, todavia falta la introduccion al metodo de elemento finito, 
sus fundamentos, se explicarán en los antecedentes, pues la teoría detrás del método 
remonta a conocimientos de algebra lineal y cálculo variacional. 
En éste primer apartado se utilizarán los conceptos básicos necesarios para una mayor 
comprensión del trabajo a elaborar, así como algunos ejemplos prácticos. 
 El segundo apartado, tratará un poco de la historia y desarrollo de los métodos, como 
fueron propuestos y las ideas detrás de ellos, es decir los antecedentes. 
En el apartado tres, se entrará de lleno a la solución de ecuaciones diferenciales parciales 
típicas, con esto se pretende mostrar la relevancia, de las EDP en el campo de la 
ingeniería por su gran aplicación. 
13 
 
Como cuarto apartado se presentarán los resultados analíticos, que serán comparados 
con los cálculos numéricos realizados en MatLab, para los diferentes tipos de ecuaciones 
diferenciales. 
En la última parte se presentarán las gráficas que arroja el MatLab, así como los códigos 
utilizados y como bonus la creación del mallado mesh. Se finalizará con una conclusión 
final sobre los alcances de los métodos analíticos y numéricos además del apéndice 
donde se encontrarán los códigos completos y las referencias citadas. 
2. Antecedentes 
 
Antes de profundizar en éste trabajo (El desarrollo de las EDP) conviene retomar el 
contexto que rodea a algunos de los temas más importantes tratados en la introducción, 
es decir sus raíces. Para ello se requirió el uso de bibliografía, misma que proporcionó 
diligentemente, el Ing. Leonardo Damián Soria Rodríguez asesor de ésta tesis y la 
profesora Nevena Perovic de la universidad técnica de München, se trata del material de 
matemáticas aplicadas II (Ecuaciones diferenciales parciales) y la de ecuaciones 
diferenciales parciales por métodos algorítmicos respectivamente. 
Ahora bien, primero se tratarán las bases para los métodos analíticos. 
2.1 Funciones ortogonales 
Es un análisis de gran utilidad, plantea que dos funciones f y g de un cierto espacio 
son ortogonales si su producto escalar 〈𝑓, 𝑔〉 es nulo, al igual que el producto punto entre 
vectores. 
La dependencia de dos funciones ortogonales depende de cómo se haya definido su 
producto escalar, es decir, de que el conjunto de funciones haya sido dotado de estructura 
de espacio. Una definición muy común de producto escalar entre funciones es: 
〈𝑓, 𝑔〉 = ∫ 𝑓(𝑥)𝑔(𝑥)𝑤(𝑥)𝑑𝑥
𝑏
𝑎
 
Con límites de integración apropiados y donde f(x) puede denotar complejo 
conjugado y w(x) es una función peso (en muchas aplicaciones se toma w(x) = 1 
Las soluciones de un problema de Sturm Liouville, es decir, las soluciones de ecuaciones 
diferenciales lineales con condiciones de borde adecuadas, pueden escribirse como una 
suma ponderada de funciones ortogonales (conocidas también como funciones propias) 
 
 
 
 
https://es.wikipedia.org/wiki/Funci%C3%B3n_matem%C3%A1tica
https://es.wikipedia.org/wiki/Producto_escalar
https://es.wikipedia.org/wiki/Espacio_prehilbertiano
https://es.wikipedia.org/wiki/Integral
https://es.wikipedia.org/wiki/Problema_de_Sturm-Liouville
14 
 
Ejemplo breve de ortogonalidad: 
f1= x2, f2= x3 para un intervalo [−1,1] y con un factor peso w(x)=1, determine si son 
funciones ortogonales. 
∫ (1) ∗ (x2) ∗ (x5)dx = 
x6
6
 valaudo de − 1 a 1 = 0
1
−1
 
 
2.2 Series de Fourier 
 
Se trata de una expansión de funciones f(x) periódicas en términos de una infinita suma 
de senos y cosenos, usando las relaciones de ortogonalidad. Las series computacionales 
de Fourier son conocidas, como análisis armónico, que es extremadamente útil para un 
determinado número de funciones arbitrarias. (González, 1997) Algunos ejemplos de 
aproximaciones sucesivas a funciones comunes usando Fourier son: 
 
Figura 1 - Series de Fourier (Wolfram Mathematica, 2016) 
Se puede apreciar que se trata de un método de ajuste de funciones, tan preciso que el 
90% de los casos, da una aproximación cercana que para una n, que tienda a infinito la 
solución, tendrá un error muy cercano al 0. 
Las series de Fourier de senos y cosenos f(x), para un intervalo de [−π, π] tiene la forma: 
𝑓(𝑥) = 
𝑎0
2
+ ∑{𝑎𝑛 cos(𝑛𝑥) + 𝑏𝑛𝑠𝑒𝑛(𝑛𝑥)}
∞
𝑛=1
 
Donde a0, an y bn son integrales de Newton o coeficientes de Fourier. Es importante notar 
aquí lo siguiente, por simetría: 
∫ sen(mx)sen(nx)dx = {
0 si n ≠ m
≠ 0 si n = m
π
−π
 
15 
 
∫ cos(mx)cos(nx)dx = {
0 si n ≠ m
≠ 0 si n = m
π
−π
 
Y donde los coeficientes de Fourier pueden ser calculados como: 
𝑎0 =
1
𝜋
∫ 𝑓(𝑥)𝑑𝑥 
𝜋
−𝜋
 
𝑎𝑛 =
1
𝜋
∫ 𝑓(𝑥) cos(𝑛𝑥) 𝑑𝑥 
𝜋
−𝜋
 
𝑏𝑛 =
1
𝜋
∫ 𝑓(𝑥) sen(𝑛𝑥) 𝑑𝑥 
𝜋
−𝜋
 
Para casos donde el intervalo sea truncado, es decir para intervalos diferentes o de 
periodo arbitrario, la estructura general es la misma, lo único que cambia es el cálculo de 
los coeficientes, que toman ahora la siguiente formulación: 
𝑎0 =
1
𝑇
∫ 𝑓(𝑥)𝑑𝑥 
𝑇/2
−𝑇/2
 
𝑎𝑛 =
2
𝑇
∫ 𝑓(𝑥) cos (
2nπ
T
x) 𝑑𝑥 
𝑇/2
−𝑇/2
 
𝑏𝑛 =
2
𝑇
∫ 𝑓(𝑥) sen (
2nπ
T
x) 𝑑𝑥 
𝑇/2
−𝑇/2
 
El caso complejo también es posible, pero no se tocará en éste trabajo, ya que las 
ecuaciones diferenciales a desarrollar no generan soluciones imaginarias. Por último pero 
no menos importante, y por facilidad de cálculos, se tiene una regla para paridad e 
imparidad de funciones: 
f(x) es par si f(x)= f(-x), o sea x2, cos(x) sólo quedará entonces el termino coseno 
f(x) es par si f(x)= -f(-x), o sea x3, sen(x) sólo quedará entonces el termino seno 
2.3 Teorema de Sturm Liouville 
 
El teorema presentado a continuación,es uno de los principales medios para la solución 
de EDP, calculadas de manera analítica, por el método de variables separables, que 
también tiene como antecedente, las series de Fourier, las cuales serán explicadas de 
manera breve en la sección de antecedentes. Sin embargo se decidió, expresar lo 
correspondiente a éste teorema en la introducción debido a su gran importancia en el 
campo, para la obtención de los valores propios (Eigen values). (R.S, 2006) 
16 
 
Las ecuaciones del tipo Sturm-Liouville, que toman su nombre de Jacques Charles 
François Sturm y Joseph Liouville, son ecuaciones diferenciales lineales de segundo 
orden de la forma: 
−
d
dx
[p(x)
dy
dx
] + q(x)y = λw(x)y 
donde p(x) y w(x) son funciones positivas semidefinidas y q(x) es real. Para el caso más 
simple, estas funciones son continuas en un intervalo finito cerrado [𝑎, 𝑏], en el que, por lo 
general, se definen unas condiciones de contorno o frontera, es decir, se concretan unos 
valores específicos que adoptan las funciones ‘’y’’. La función w(x) es llamada función 
peso. 
Por el otro lado el valor de λ, no se especifica, pues se desea encontrar los valores de 
ella, para los que exista una solución no trivial de la ecuación que satisfaga las 
condiciones de frontera se denomina el problema (S-L). 
Esas lambdas son conocidas como los valores propios del problema de S-L que se 
plantea en la ecuación de la hoja anterior, en conjunto con las condiciones de frontera. 
Las soluciones correspondientes son las funciones propias o las eigenfunciones o los 
eigenvectores del problema. Éste tema, como se había mencionado anteriormente, es de 
suma importancia en matemáticas aplicadas, donde los problemas S-L ocurren muy 
comúnmente, particularmente al resolver ecuaciones diferenciales parciales con 
separación de variables. 
Ejemplo: 
Se tiene una ecuación diferencial de segundo orden, definida por 
𝒚’’ + 𝝀𝒚 = 𝟎 
Aplicando los conociminetos previos de ecuaciones diferenciales de orden tenemos los 
casos posibles: λ<0, λ=0 y λ>0. 
Para el caso λ<0, la solución quedaria descrita como: 
𝑦(𝑥) = 𝐶1
√−λ2𝑥 + 𝐶2
−√−λ2𝑥 
Para el caso λ=0, la solución quedaria descrita como: 
𝑦(𝑥) = 𝐶1 + 𝐶2𝑥 
Y por último para λ>0, la solución quedaria descrita como: 
𝑦(𝑥) = 𝐶1𝑐𝑜𝑠(√λ𝑥) + 𝐶2𝑠𝑒𝑛(√λ𝑥) 
Entonces, puede notarse para el caso en que lambda es igual a 0, que el unico camino 
viable es 𝐶1 = 𝐶2 = 0, es decir la solución trivial. En caso contrario, a lambdas menores y 
mayores que cero, no hay valor propio que las haga cero. 
https://es.wikipedia.org/wiki/Funci%C3%B3n_real
https://es.wikipedia.org/wiki/Trivial_(matem%C3%A1tica)
17 
 
Por lo tanto éstas son las que generan las funciones propias, mismas que contienen de 
manera explicita las series de Fourier. 
Si λ>0 
𝑦(0) = 𝐶1𝑐𝑜𝑠(√λ(0)) + 𝐶2𝑠𝑒𝑛(√λ(0)) 
𝑦(0) = 𝐶1𝑐𝑜𝑠(0) + 𝐶2𝑠𝑒𝑛(0) 
Por lo tanto 𝐶1 = 0 , pero si 
𝑦(1) = 𝐶2𝑠𝑒𝑛(√λ) = 0 
𝑠𝑒𝑛(√λ) = 0 que da λn = nπ 
Función propia: y = Cnsen(√nπx) 
2.4 Ecuaciones de Navier-Stokes 
Se habla sin duda, de una de las ecuaciones que revolucionó al mundo, debido a su 
aplicación en la reología, o mecánica de fluidos. Éstas ecuaciones toman su nombre de 
dos investigadores: Claude-Louis Navier y George Gabriel Stokes, se trata de un 
conjunto de ecuaciones en derivadas parciales no lineales que describen el movimiento 
de un fluido, sus rotaciones y mecánica. Son hoy en día aplicables a fenómenos 
atmosféricos, corrientes oceánicas y el flujo alrededor de vehículos o proyectiles y, en 
general a cualquier fenómeno en el que se involucren fluidos newtonianos. 
Fueron obtenidas, aplicando los principios de conservación de la mecánica y 
la termodinámica a un volumen fluido. Haciendo esto se obtiene la llamada formulación 
integral de las ecuaciones. Para su forma diferencial, se manipulan aplicando ciertas 
consideraciones, principalmente aquellas en las que los esfuerzos tangenciales guardan 
una relación lineal con el gradiente de velocidad, la famosa ley de viscosidad de Newton, 
obteniendo de esta manera, la formulación diferencial que generalmente es más útil para 
la resolución de los problemas que se plantean en la mecánica de fluidos. 
No se dispone de una solución general para este conjunto de ecuaciones, y salvo ciertos 
tipos de flujo y situaciones muy concretas, no es posible hallar una solución analítica; por 
lo que en muchas ocasiones es preciso recurrir al análisis numérico para determinar una 
solución aproximada. 
La ecuación tiene en el particular, los siguientes tres casos, para viscosidad despreciable, 
es decir µ=0: 
𝜌 (
𝜕𝑣𝑥
𝜕𝑡
+ 𝑣𝑥
𝜕𝑣𝑥
𝜕𝑥
+ 𝑣𝑦
𝜕𝑣𝑥
𝜕𝑦
+ 𝑣𝑧
𝜕𝑣𝑥
𝜕𝑧
) = µ [
𝜕2𝑣𝑥
𝜕𝑥2
+
𝜕2𝑣𝑥
𝜕𝑦2
+
𝜕2𝑣𝑥
𝜕𝑧2
] −
𝜕𝑃
𝜕𝑥
− 𝜌𝑔𝑥 
𝜌 (
𝜕𝑣𝑦
𝜕𝑡
+ 𝑣𝑥
𝜕𝑣𝑦
𝜕𝑥
+ 𝑣𝑦
𝜕𝑣𝑦
𝜕𝑦
+ 𝑣𝑧
𝜕𝑣𝑦
𝜕𝑧
) = µ [
𝜕2𝑣𝑦
𝜕𝑥2
+
𝜕2𝑣𝑦
𝜕𝑦2
+
𝜕2𝑣𝑦
𝜕𝑧2
] −
𝜕𝑃
𝜕𝑦
− 𝜌𝑔𝑦 
https://es.wikipedia.org/wiki/Fluido
18 
 
𝜌 (
𝜕𝑣𝑧
𝜕𝑡
+ 𝑣𝑥
𝜕𝑣𝑧
𝜕𝑥
+ 𝑣𝑦
𝜕𝑣𝑧
𝜕𝑦
+ 𝑣𝑧
𝜕𝑣𝑧
𝜕𝑧
) = µ [
𝜕2𝑣𝑧
𝜕𝑥2
+
𝜕2𝑣𝑧
𝜕𝑦2
+
𝜕2𝑣𝑧
𝜕𝑧2
] −
𝜕𝑃
𝜕𝑧
− 𝜌𝑔𝑧 
y la ecuación de continuidad adquiere la forma siguiente: 
𝜕𝑣𝑥
𝜕𝑥
+
𝜕𝑣𝑦
𝜕𝑦
+
𝜕𝑣𝑧
𝜕𝑧
= 0 
Como se mencionó anteriormente, la ecuación no tiene una solución única, pues depende 
en gran manera del problema que se esté abordando, y las condiciones a la frontera que 
se le apliquen. Un caso particular y derivado de éste modelo, es el caso de la ecuación de 
difusión, que en su estado estacionario es la de Laplace. 
 Dicha ecuación ha sido estudiada anteriormente en problemas de transferencia de masa 
y es muy socorrida en muchos problemas de fenómenos de transporte y polímeros. 
 
2.5 Antecedentes - Método de elemento finito (MEF ó FEM) 
 
Para la descripción de los fenómenos que ocurren en la naturaleza, se presentan muchas 
veces modelos en forma de ecuaciones diferenciales parciales. Algunas de las cuales no 
pueden ser resueltas por los métodos conocidos, sin embargo, es de suma importancia 
tener una opción, que permita resolver dichos problemas. 
Uno de esas alternativas, que se presenta cada vez con mayor frecuencia, es el método 
del elemento finito (MEF en español) que es conocido por sus buenos resultados y 
aplicaciones prácticas, no sólo en el área de la ingeniería mecánica, sino también en la 
Ing. Química, como se propone como eje central de ésta tesis. El método consiste 
principalmente en transformar una ecuación diferencial, en un sistema de ecuaciones 
algebraicas a partir de una discretización, dicho de otro modo, separarla en varias 
particiones o elementos finitos del dominio de definición de la función, para luego, usando 
algún método aproximado, obtener la solución en cada uno de los elementos, es decir se 
tomará la aportación de cada uno de ellos y se sumará, con el propósito de obtener un 
resultado con un margen de error muy cercano al cero. 
Éste tema será expuesto más adelante (Último método de solución para EDP) pues antes 
se requiere de algunos conceptos y trabajos como: 
2.6 Independencia lineal 
Es un tema de gran importancia para el algebra lineal, trata basicamente de conjuntos de 
vectores, donde si son linealmente independientes, un vector no pueda ser escrito como 
combinación lineal de otro (Ya sea que sumen dos y de un tercero o un multiplo de 
alguno) Un ejemplo en R3, el conjunto o base canónica e3 es decir (1,0,0), (0,1,0), (0,0,1) 
son linealmente independientes unos de otros, caso contrario pasa cuando se tiene (4,-
2,2), (1,0,1), (5,-2,3) pues el tercer elemento es la suma de los primeros dos, lo cual 
indica que no es una base, ya que uno de sus vectores es linealmente dependiente. 
19 
 
Una manera para comprobar si los vectores dentro de una matriz son todos lienalmente 
indepenientes, es haciendo submatrices, y entonces observar si el determinante calculado 
es igual a 0; entonces se trata de un sistema deecuaciones con un grado de libertad, por 
ende el sistema es inestable ya que se tienen vectores como resultado de combinaciones 
lineales. 
La condicion que se debe satisfacer es: 
𝛼1𝑉1 + 𝛼2𝑉2 + 𝛼3𝑉3 + ⋯ + 𝛼𝑛𝑉𝑛 = 0 α ϵ R, donde las alfas 𝛼1 = 𝛼2 = 𝛼3 = 𝛼𝑛 = 0 
2.7 Métodos variacionales de Rayleight Ritz 
 
Método implementado originalmente, para las ecuaciones del tipo Lagrange y Euler, 
asociadas con una función Lagrangiana. La técnica se basa en representar los 
movimientos de un sistema dinámico, mediante la combinación lineal de una serie de 
perturbaciones linealmente independientes, es decir a las fuerzas aplicadas sobre un 
cuerpo o sustancia que cumplan las condiciones de contorno geométricas: 
𝑦(𝑥, 𝑡) = ∑ 𝑁𝑗(𝑥)𝑣𝑗(𝑡)
𝑛
𝑗−1
 
Donde 𝑁(𝑥) son las perturbaciones o deformaciones, o funciones de forma de su nombre 
en ingles (N- Shape functions), y las 𝑣𝑗son los coeficientes escalares asociados a cada 
una de las deformaciones. La condición que debe ser mantenida es que las N(x) sean 
linealmente independientes unas de otras, no implica que alguna de ellas no pueda ser 
representada como combinacion lineal de las otras. (J., 2002) 
De hecho el método de elemento finito, es un caso particular de éste problema de 
Rayleight, en el que las funciones de forma 𝑁(𝑥) o shape functions, tienen las 
propiedades anteriormente mencionadas, que facilitarán las matrices de masa, rigidez y la 
de vectrores fuerza o resitencia. 
Éste método a diferencia del MEF, es que Ritz permite usar cualquier combinacion de 
funciones de forma linealmente independientes y que cumpla condiciones de contorno, lo 
que permite que los coeficientes de la ecuacion no tengan que presentar un giro concreto. 
El grado de aproximación de las deformaciones, fluctuará dependiendo el número de 
términos empleados, es decir del número de elementos, o de coordenadas de cada uno 
de ellos. (J., 2002) 
La solución se determina a través de 
𝑆 = ∫ 𝑑𝑥𝐿(𝑌, 𝑌′, 𝑥)
𝑏
𝑎
 
20 
 
Donde Rayleight, sustituira la solucion exacta es decir la Y(x) por una aproximada y(x), 
proveniente de los parámetros 𝑘1, 𝑘2, 𝑘3… que son los elementos locales de una matriz, 
que se estudiará más adelante, ahora S será: 
𝑑𝑆[𝑦(𝑥, 𝑘1, 𝑘2, 𝑘3,… ] = 0 
O sea 
𝜕𝑆
𝜕𝑘1
= 
𝜕𝑆
𝜕𝑘2
= ⋯ = 0 se planteará un sistema de n ecuaciones, con n incognitas, 
para obtener los parámetros k. 
2.8 Mallado 
Funciones del mallado: 
1. Transformar lo continuo en segmentos o elementos, mediante líneas o superficies 
imaginarias en un número de elementos finitos. (iberisa, 2002) 
2. Los elementos están conectados entre sí, mediante un número discreto de puntos o 
“nodos”, situados en sus contornos. (iberisa, 2002) 
3. Se toma un conjunto de funciones que definen de manera única, el campo de 
desplazamientos dentro de cada “elemento finito” en función de los desplazamientos 
nodales de dicho elemento. Por ejemplo el campo de desplazamientos dentro de un 
elemento lineal de dos nodos podría venir definido por: U = N1U1 + N2U2, siendo N1 y N2 
las funciones comentadas (funciones de forma) y u1 y u2 los desplazamientos en el nodo 
1 y en el nodo 2. (iberisa, 2002) 
4. Las funciones de desplazamiento o soluciones U, definen el estado de deformación del 
elemento en función de los desplazamientos nodales. 
5. Nos ayuda a determinar cuáles son las fuerzas concentradas en los nodos, resultando así 
una relación entre fuerzas y desplazamientos de la forma F = K·u, que como vemos es 
similar a la del cálculo matricial. (iberisa, 2002) 
Algunas reglas para el mallado que se deben seguir: 
 Para el análisis de tensiones se requiere una malla más densa. estudios 
térmicos requieren mallas aún más densas. (iberisa, 2002) 
 Los nodos tienen que tener una localización estratégica en general donde sea 
necesario obtener información de desplazamientos, tensiones o temperaturas. 
 Recordar que los nodos son la única fuente de información (la información que 
dan los elementos se deriva de los nodos). (iberisa, 2002) 
21 
 
 
Figura 2 - Nodos 
 Si es posible, usar mallado uniforme (Igual separación entre nodos). En regiones de 
transición de alta a baja densidad de malla, no cambiar las dimensiones de los 
elementos adyacentes por un factor mayor de 2. (iberisa, 2002) 
 
Figura 3 - Mallado 
 Mallar con cuadriláteros, en lugar de elementos triangulares. Usar únicamente 
elementos triangulares sólo para transiciones de malla o por exigencias geométricas. 
(iberisa, 2002) 
 
Figura 4 – Transición 
 
 
 
 
22 
 
3 Objetivos 
 
3.1 Objetivos generales 
 
 Desarrollar una metodología para resolver ecuaciones diferenciales parciales en 
ecuaciones de amplio espectro (Momentum, energía y masa) mediante la técnica analítica 
de variables separables y el método numérico de elemento finito y diferencias finitas de 
Euler. 
 Desarrollar de manera detallada los algoritmos o pasos a seguir para resolver las EDP, en 
una, dos y tres dimensiones. 
 Implementar en un programa computacional la solución numérica del método de elemento 
finito. 
 Demostrar las ventajas y desventajas de ambos modelos, así como una comparativa 
entre los mismos, con el fin de determinar cuál de ambos aplica mejor, bajo determinadas 
circunstancias. 
 
3.2 Objetivos particulares 
 
 Describir el método analítico y numérico detalladamente. 
 Tratar las ecuaciones diferenciales de manera algorítmica en sus tres diferentes tipos: 
Elípticas, hiperbólicas y parabólicas. 
 Analizar las diferencias entre las ecuaciones diferenciales parciales de una, dos y 3 
dimensiones. 
 Calcular la solución numérica, con la paquetería ofrecida por el software de MatLab, es 
decir gráficas, cálculos matriciales y ciclos. 
 Expresar de manera explícita, mediante ejemplos diversos, la solución numérica y su 
implementación computacional para resolver ecuaciones parciales en una, dos y 3 
dimensiones. 
 Demostrar que la aplicación del método numérico, facilita en gran escala los cálculos para 
la solución de ecuaciones diferenciales parciales en cualquier dimensión que se trate. 
 
Para alcanzar los objetivos deseados, se desarrollarán en los capítulos 2, 3 y 4 los 
métodos numéricos y analíticos, así como los algoritmos para la solución de ecuaciones 
diferenciales parciales. En el capítulo último se mostrarán las aplicaciones físicas e 
ingenieriles del trabajo realizado, así como sus alcances. 
 
 
 
23 
 
4 Desarrollo del método analítico 
 
4.1 Caso I – Perfil de velocidades 
 
 ∇𝟐𝐕 =
𝛛𝟐𝐕
𝐝𝐱𝟐
 + 
𝛛𝟐𝐕
𝐝𝐲𝟐
= 0 …. (1) 
Para las condiciones a la frontera dadas: 
𝑉(0, 𝑦) = 𝑔(𝑦) … (2) y 𝑉(𝐿, 𝑦) = 0 … (4) 
 y para 𝑉(𝑥, 0) = 0 … (3) con 𝑉(𝑥, 𝐻) = 0 … (5) 
Se propone una solución por variables separables: 
Proponer una función basada en principio de superposición: 
𝑉(𝑥, 𝑦) = 𝛷(𝑥) ∗ 𝛳(𝑦) … (6) 
Se sustituye en la ecuación original, o sea la 6: 
∂2Φ(x)∗ϴ(y)
dx2
+
∂2Φ(x)∗ϴ(y)
dy2
= 0… (7) 
𝛷’’(𝑥) ∗ 𝛳(𝑦) + 𝛷(𝑥) ∗ 𝛳’’(𝑦) = 0 … (8) 
Se divide entre la ecuación 6 y se iguala a una constante de separación: 
−
Φ′′(x)
Φ(x)
=
ϴ′′(y)
ϴ(y)
= −λ2 … (9) 
Se obtiene el siguiente sistema de ecuaciones diferenciales ordinarias: 
 𝛷’’(𝑥) − 𝜆2𝛷(𝑥) = 0 … (10) 𝛷(𝐿) = 0 … (12) 
 𝛳’’(𝑦) + 𝜆2 𝛳 (𝑦) = 0 … (11) 𝛳 (0) = 0 … (13) 𝑦 𝛳 (𝐻) = 0 … (14) 
Se resuelve cada una de las ecuaciones con las condiciones a la frontera 13 y 14, éstas 
son ecuaciones diferenciales ordinarias de segundo grado, se puede optar por reducción 
de orden o por el método de la exponencial y senos y cosenos para las raíces complejas. 
𝛳 (𝑦) = 𝐴𝑠𝑒𝑛(𝜆𝑦) + 𝐵𝑐𝑜𝑠(𝜆𝑦) 
𝛳 (0) = 𝐵 = 0 
𝛳 (𝐻) = 𝐴𝑠𝑒𝑛(𝐻𝑦) 
Se expresa entonces el valor propio como como: 
𝐻𝜆 = 𝑛𝜋 por lo tanto, 𝜆 = 𝑛𝜋
𝐻
 … (15) 
24 
 
𝛳 (𝑦) = 𝐴𝑛𝑠𝑒𝑛 (
𝑛𝜋
𝐻
𝑦) … (16) 
Solución particularde la coordenada y, se sigue la misma metodología para la ecuación 
diferencial ordinaria de la coordenada x: 
𝛷(𝑥) = 𝐶1𝑐𝑜𝑠ℎ (
𝑛𝜋
𝐻
𝑥) + 𝐶2𝑠𝑒𝑛ℎ (
𝑛𝜋
𝐻
𝑥) … (17) 
𝛷(𝐻) = 0 
𝛷(𝑥) = 𝐶1𝑐𝑜𝑠ℎ (
𝑛𝜋
𝐻
∗ 𝐿) + 𝐶2𝑠𝑒𝑛ℎ (
𝑛𝜋
𝐻
∗ 𝐿) … (18) 
Lo cual no sirve, pero 17 se puede expresar utilizando 18 de la siguiente manera: 
𝛷(𝑥) = 𝐶1𝑐𝑜𝑠ℎ (
𝑛 ∗ 𝜋 ∗ (𝑥 − 𝐿)
𝐻
) + 𝐶2𝑠𝑒𝑛ℎ (
𝑛 ∗ 𝜋 ∗ (𝑥 − 𝐿)
𝐻
) … (19) 
𝛷(𝐿) = 𝐶1𝑐𝑜𝑠ℎ (
𝑛 ∗ 𝜋 ∗ (𝐿 − 𝐿)
𝐻
) + 𝐶2𝑠𝑒𝑛ℎ (
𝑛 ∗ 𝜋 ∗ (𝐿 − 𝐿)
𝐻
) … (20) 
Usando 𝛷(𝑥) = 𝐶1 = 0, entonces 
𝛷(𝑥) = 𝐶2𝑠𝑒𝑛ℎ (
𝑛 ∗ 𝜋 ∗ (𝑥 − 𝐿)
𝐻
) … (21) 
Se puede calcular entonces la ecuación para la coordenada x. 
Una vez con ambas soluciones calculadas de manera independiente se sustituyen en 6, 
para obtener la solución general. 
𝑉(𝑥, 𝑦) = 𝐴𝑛 ∗ 𝑠𝑒𝑛ℎ (
𝑛𝜋
𝐻
∗ 𝑦) 𝐶2𝑠𝑒𝑛ℎ (
𝑛 ∗ 𝜋 ∗ (𝑥 − 𝐿)
𝐻
)… (22) 
Ese es el perfil de velocidades, en su forma general, expresada en términos de x y y. 
Mediante del principio de superposición, se calculan las constantes enésimas, dado que 
se trata de una multiplicación, se puede calcular, se puede reescribir la ecuación como 
una serie, que va desde n igual a uno hasta infinito, y donde la Cn, puede representar una 
nueva constante, que es a su vez la multiplicación de An por C2. 
𝑉(𝑥, 𝑦) = ∑ 𝐶𝑛𝑠𝑒𝑛 (
𝑛𝜋
𝐻
∗ 𝑦) ∗ 𝑠𝑒𝑛ℎ (
𝑛 ∗ 𝜋 ∗ (𝑥 − 𝐿)
𝐻
)
∞
𝑛=1
 … (23) 
Aplicando la segunda condición a la frontera se obtiene una nueva función de y: 
𝑔(𝑦) = ∑ 𝐶𝑛 ∗ 𝑠𝑒𝑛 (
𝑛𝜋
𝐻
∗ 𝑦) ∗ 𝑠𝑒𝑛ℎ (
𝑛 ∗ 𝜋 ∗ (𝑥 − 𝐿)
𝐻
)
∞
𝑛=1
 … (24) 
25 
 
Se procede a usar el principio de ortogonalidad, para el cálculo de las constantes (revisar 
antecedentes) y series de Fourier para entender la parte conceptual: 
 ∫ 𝑔(𝑦) ∗ 𝑠𝑒𝑛(
𝑛 ∗ 𝜋
𝐻
∗ 𝑦)
𝐻
0
𝑑𝑦 = 𝐶𝑛 ∗ ∫ 𝑠𝑒𝑛
2 (
𝑛 ∗ 𝜋
𝐻
∗ 𝑦)
𝐻
0
 𝑠𝑒𝑛ℎ(
−𝑛 ∗ 𝜋 ∗ 𝐿
𝐻
)𝑑𝑦 … (25) 
Para encontrar la constante Cn, se propone una solución lineal y una trigonométrica: 
𝐶𝑛 = 
∫ 𝑔(𝑦) ∗ 𝑠𝑒𝑛(
𝑛 ∗ 𝜋
𝐻 ∗ 𝑦)
𝐻
0
𝑑𝑦 
𝑠𝑒𝑛ℎ (
−𝑛 ∗ 𝜋 ∗ 𝐿
𝐻 ) ∗ ∫ 𝑠𝑒𝑛2 (
𝑛 ∗ 𝜋
𝐻 ∗ 𝑦)
𝐻
0
 𝑑𝑦
 … (26) 
𝑉(𝑥, 𝑦) = ∑
∫ g(y)∗sen(
n∗π
H
∗y)
H
0 dy
senh(
−n∗π∗L
H
)∗∫ sen2(
n∗π
H
∗y)
H
0 dy
∞
𝑛=1 … (27) 
4.2 Caso II – Ecuación de transferencia de momentum, con cambio en las C.F 
 
Para la ecuación diferencial parcial, se trata de un caso similar al anterior, sólo que esta 
vez, cambian las condiciones a la frontera para la coordenada x, generando una solución 
muy similar, sin embargo diferente a la calculada anteriormente: 
 𝛁𝟐𝑽 =
𝝏𝟐𝑻
𝝏𝒙𝟐
+
𝝏𝟐𝑻
𝝏𝒚𝟐
= 𝟎 … (28) 
Para las condiciones a la frontera de: 
𝑉(0, 𝑦) = 0 … (29) ; 𝑉(𝐿, 𝑦) = 0 … (30) 
 𝑉(𝑥, 0) = 0 … (31) ; 𝑉(𝑥, 𝐻) = 𝑓(𝑥) … (32) 
Nuevamente, se propone una solución por variables separables, para ello se ocupará 
nuevamente el principio de superposición: 
𝑉(𝑥, 𝑦) = 𝜁(𝑥) ∗ 𝜑(𝑦) …(33) 
Se sustituye nuevamente en la ecuación original 
∂2𝜁(x)∗𝜑(y)
dx2
+
∂2𝜁(x)𝜑(y)
dy2
= 0 … (34) 
𝜁’’(x) 𝜑 (y) + +𝜁(x)𝜑’’(y) = 0 … (35) 
Se divide entre la ecuación 6 y se iguala a una constante de separación: 
−
𝜁′′(x)
𝜁(x)
=
𝜑′′(y)
𝜑(y)
= −𝜆2…(36) 
Se obtiene el siguiente sistema de ecuaciones diferenciales ordinarias: 
 𝜁"(x) − 𝜆2 𝜁 (x) = 0… (37) 𝜑 (𝐿) = 0 … (39) 
26 
 
 𝜑’’(𝑦) + 𝜆2 𝜑(𝑦) = 0 … (38) 𝜁 (0) = 0 … (40) 𝑦 𝜁 (𝐻) = 0 … (41) 
Nuevamente se resuelven las ecuaciones con las condiciones a la frontera 40 y 41, que 
generarán, dos nuevas ecuaciones diferenciales ordinarias, se propone el método de 
senos y cosenos para las raíces complejas. 
 ϴ (y) = Asen(λx) + Bcos(λx) …(42) 
ϴ (0)= B =0 
ϴ (H)= Asen(Hy) 
Se expresa entonces el valor propio como: 
Hλ = nπ por lo tanto, 𝜆 = 𝑛𝜋
𝐻
 … (42) 
ϴ (y) = Ansen (
nπ
H
x)… (43) 
Solución particular de la coordenada y, con la misma metodología que se usó para la 
ecuación diferencial ordinaria de la coordenada x, explicada en el caso anterior: 
Φ(x) = 𝐶3cosh (
nπ
H
y) + 𝐶4senh (
nπ
H
𝑦) … (44) 
𝛷(𝐻) = 0… (44𝑎) 
Φ(x) = 𝐶3cosh (
nπ
H
𝐿) + 𝐶4senh (
nπ
H
𝐿) … (45) 
Lo cual no sirve, pero 44 se puede expresar utilizando 45 de la siguiente manera: 
𝛷(𝑥) = 𝐶3𝑐𝑜𝑠ℎ (
𝑛𝜋(𝑦−𝐿)
𝐻
) + 𝐶4𝑠𝑒𝑛ℎ (
𝑛𝜋(𝑦−𝐿)
𝐻
) … (46) 
𝛷(𝐿) = 𝐶3𝑐𝑜𝑠ℎ (
𝑛𝜋(𝐿−𝐿)
𝐻
) + 𝐶4𝑠𝑒𝑛ℎ (
𝑛𝜋(𝐿−𝐿)
𝐻
) …(47) 
Usando 𝛷(𝑥) = 𝐶3 = 0, entonces 
Φ(x) = 𝐶4senh (
nπ(y−L)
H
)… (48) 
Se puede calcular entonces la ecuación para la coordenada x. 
Una vez con ambas soluciones calculadas de manera independiente se sustituyen en 33, 
para obtener la solución general. 
V(x, y) = Bnsenh (
nπ
H
x) 𝐶3senh (
nπ(y−L)
H
)… (49) 
Ese es el perfil de velocidades, en su forma general, expresada en términos de x y y. 
27 
 
A través del principio de superposición, se obtienen las constantes enésimas como en el 
caso anterior, dado que se trata de una multiplicación, se reescribe la ecuación como una 
serie, que va desde n igual a uno hasta infinito, y donde la Dn, puede representar una 
nueva constante, que es a su vez la multiplicación de Bn por 𝐶4. 
V(x, y) = ∑ Dnsen (
nπ
H
x) senh (
nπ(y−L)
H
)∞𝑛=1 …(50) 
Aplicando la segunda condición a la frontera se obtiene una nueva función de y: 
f(x) = ∑ Dnsen (
nπ
H
x) senh (
nπ(y−L)
H
)∞𝑛=1 …(51) 
Se procede a usar el principio de ortogonalidad, para el cálculo de las constantes (revisar 
antecedentes) y series de Fourier para entender la parte conceptual: 
∫ f(x) ∗ sen(
n∗π
H
x)
H
0
𝑑𝑥 = Dn∫ sen2 (nπ
H
x)
H
0
senh(
−nπL
H
)dx … (52) 
Para encontrar la constante Cn, proponer solución lineal y trigonométrica 
𝐷𝑛 = 
∫ 𝑓(𝑥)𝑠𝑒𝑛(
𝑛𝜋
𝐻 𝑥)
𝐻
0
𝑑𝑥 
𝑠𝑒𝑛ℎ (
−𝑛𝜋𝐿
𝐻 ) ∗ ∫ 𝑠𝑒𝑛
2 (
𝑛𝜋
𝐻 𝑥)
𝐻
0
 𝑑𝑥
 … (53) 
𝑉(𝑥, 𝑦) = (∑
∫ f(x)∗sen(
nπ
H
x)
H
0
dx
senh(
−nπL
H
)∫ 𝑠𝑒𝑛2(
nπ
H
x)
H
0
 dx
∞
𝑛=1 ) [sen (
nπ
H
x) senh (
nπ(y−L)
H
)] …(54) 
4.3 Caso III – Ecuación de difusión de energía para el caso de 2 dimensiones 
 
La ecuación de difusión, muy socorrida en la transferencia de energía, se puede obtener 
de las ecuaciones del balance de energía, desarrollada en la parte de antecedentes, tiene 
la forma: 
𝛛𝐓
𝛛𝐭
= −𝛂𝟐𝛁𝟐𝐓 
∂T
∂t
= (
∂2T
∂x2
+
∂2𝑇
∂y2
 ) 𝛼2 … (55) 
Con condiciones a la frontera: 
𝑇(0, 𝑦, 𝑡) = 0 … (56) ; 𝑇(𝑎, 𝑦, 𝑡) = 0 … (57) ; 𝑇(𝑥, 𝑦, 0) = 𝜔(𝑥, 𝑦) … (58) 
𝑇(𝑥, 0, 𝑡) = 0 … (59) ; 𝑇(𝑥, 𝑏, 𝑡) = 0 … (60) 
Se usa como es tradición el método de variables separables, es decir se propone una solución 
para el caso de 3 variables: 
𝑇(𝑥, 𝑦) = 𝛷(𝑥) ∗ 𝛳(𝑦) ∗ 𝜗(𝑡) … (61) 
28 
 
Igual que los dos casos se iguala a una constante de separación: 
 𝜗′(t)
𝜗
 = (
Φ′′(x)
Φ(x)
+
ϴ′′(y)
ϴ(y)
 ) 𝛼2 = −λ 2 …(62) 
Se considera que la -λ2 puede ser dividida sobre el 𝛼2 para comprimirlas en una nueva 
variable −𝛽2 negativo, que conservará el signo de la lambda. 
𝜗′(t)
𝜗
= (
Φ′′(x)
Φ(x)
+
ϴ′′(y)
ϴ(y)
 ) =
−λ 2
𝛼2
= −𝛽2 …(63) 
Y ésta sustituida como en los casos anteriores, genera un nuevo modelo de sistema de 
ecuaciones diferenciales ordinarias de segundo grado, sólo que ésta vez, se originan 3 
ecuaciones, para las variables t, x y y: 
𝜗’(𝑡) + 𝛽2𝑡 = 0 … (64) 
 𝛷′′(𝑥) + 𝛾2 𝛷 = 0 … (65) 
𝛳′′(𝑦) + 𝛿2𝛳(𝑦) = 0 … (66) 
 La solución a las mismas, se obtiene por despeje de la ecuación 64, 
𝛽2= γ2+𝛿2 …(67) 
𝜗(t)=A𝑒−𝛽2t …(68) 
Φ(x) = B cos(ρx) + Csen(ρx) …(69) 
ϴ(y) = E ∗ sen(δy) + D ∗ cos (δy) …(70) 
β2 =
n2π2
a2
+ 
𝑚2π2
b2
 …(71) 
λ2 = α2 ∗ (
n2π2
a2
+ 
m2π2
b2
) …(72) 
Se tiene una constante que depende de dos contadores n y m o sea la Fnm 
T(x, y, t) = ∑ ∑ Fnm
∞
m=1
∞
n=1 𝑒
−α2∗(
n2π2
a2
+ 
m2π2b2
)∗𝑡
 sen (
nπ
a
∗ x) sen(
mπ
b
y) …(73) 
Aplicando las condiciones a la frontera: 
T(x, y, t) = ∑ ∑ Fnm ∗ sen (
n∗π
a
∗ x)∞m=1
∞
n=1 sen(
m∗π
b
∗ y) = ω(x, y) …(74) 
 ∫ ∫ sen (nπ
a
x) sen (
mπ
b
y)
b
0
𝑤(𝑥, 𝑦)dxdy 
𝑎
0
= Fnm ∫ sen
2 (
nπ
a
x)
a
0
dx ∫ sen2 (
mπ
b
y) dy. . (75)
b
0
 
Fnm = ∑ ∑
∫ ∫ sen(
nπ
a
x)sen(
mπ
b
y)∗
b
0
𝑤(𝑥,𝑦)dxdy 
𝑎
0
 
∫ sen2(
nπ
a
x)
a
0
dx ∫ sen2(
mπ
b
y)dy
b
0
∞
𝑚=1
∞
𝑛=1 …(76) 
29 
 
T(x,y,t ) =∑ ∑
∫ ∫ sen(
nπ
a
∗x)sen(
mπ
b
y)
b
0
𝑤(𝑥,𝑦)dxdy 
𝑎
0
 
∫ sen2(
nπ
a
x)
a
0
dx∫ sen2(
mπ
b
y)dy
b
0
∞
𝑚=1
∞
𝑛=1 *𝑒
−α2(
n2π2
a2
+ 
m2π2
b2
)𝑡* sen (nπ
a
x) sen(
mπ
b
y) …(77) 
Solución particular, expresada en términos de una condición inicial general, dadas por la w(x,y) 
T(x,y,t)= ∑ ∑ 4
𝑎𝑏
∫ ∫ sen (
nπ
a
x) ∗ sen (
mπ
b
y)
b
0
𝑤(𝑥, 𝑦)dxdy ∗ 𝑒
−α2∗(
n2π2
a2
+ 
m2π2
b2
)∗𝑡
 sen (
nπ
a
∗ x) sen(
mπ
a
y) 
𝑎
0
∞
m=1
∞
n=1 …(78) 
4.4 Caso IV – Ecuación de difusión de energía para el caso de 3 dimensiones 
 
Lo mismo se aplicará con el caso de tres dimensiones, sólo que ahora, es de suponerse que la ecuación ya no involucrará dos sumatorias, 
sino más bien tres series, por tener una tercera parcial con respecto a z en ésta ocasión, la ecuación tiene la forma ya conocida: 
∂T
∂t
= −α2∇2U ∂T
∂t
= (
∂2T
∂x2
+
∂2𝑇
∂y2
+
∂2𝑇
∂z2
) 𝛼2 …(79) 
Con condiciones a la frontera previamente usadas: 
𝑇(0, 𝑦, 𝑡) = 0 … (56) ; 𝑇(𝑎, 𝑦, 𝑡) = 0 … (57) ; 𝑇(𝑥, 𝑦, 𝑧, 0) = 𝜔(𝑥, 𝑦, 𝑧) … (58) 
𝑇(𝑥, 0, 𝑡) = 0 … (59) ; 𝑇(𝑥, 𝑏, 𝑡) = 0 … (60) 
Se usa como es tradición, el método de variables separables, es decir se propone una solución para el caso de 3 variables: 
𝑇(𝑥, 𝑦) = 𝛷(𝑥) ∗ 𝛳(𝑦) ∗ 𝛹(𝑧) ∗ 𝜗(𝑡) … (61) 
Igual que los dos casos se iguala a una constante de separación: 
 𝜗′(t)
𝜗
 = (
Φ′′(x)
Φ(x)
+
ϴ′′(y)
ϴ(y)
+ 
Ψ′′(z)
Ψ(z)
)𝛼2 = −𝜆2 …(62) 
30 
 
Se considera que la -λ2 puede ser dividido sobre el 𝛼2 a manera de comprimirlas en una 
nueva variable −𝛽2 negativo pues conserva el signo de la lambda. 
𝜗′(t)
𝜗
 = (Φ′′(x)
Φ(x)
+
ϴ′′(y)
ϴ(y)
+
Ψ′′(z)
Ψ(z)
)= −𝜆2
𝛼2
 = −𝜅2 …(63) 
Y ésta sustituida como antes, lo que generá un nuevo modelo de sistema de ecuaciones 
diferenciales ordinarias de segundo grado, sólo que ésta vez, se originan 4 en lugar de 
dos, para las variables: t, x, y, z: 
𝜗’(t)+ 𝛽 2*t=0 …(64) 
 Φ′′(x)+γ2Φ=0 …(65) 
ϴ′′(y) + 𝛿2ϴ(y) =0 …(66) 
Ψ′′(z) + 𝜂2Ψ = 0 …(67) 
 La solución a las mismas, se obtiene por despeje de la ecuación 64, 
𝜅2 = 𝛾2 + 𝛿2 + 𝜂2 …(68) 
𝜗(t) = A𝑒−𝜅
2t …(69) 
Φ(x) = B cos(ρx) + Csen(ρx) …(70) 
ϴ(y) = Esen(δy) + Dcos (δy) …(71) 
Ψ(z) = Fsen(𝜂y) + Gcos (𝜂y) …(72) 
𝜅2 =
n2π2
a2
+ 
n2π2
b2
 …(73) 
λ2 = α2 (
n2π2
a2
+ 
m2π2
b2
+
l2π2
c2
 ) …(74) 
Se tiene una cte que depende de dos contadores n, m y l, o sea la Hnml 
T(x, y, z, t) = ∑ ∑ ∑ Hnml
∞
l=1
∞
m=1
∞
𝑛=1 𝑒
−α2(
n2π2
a2
+ 
m2π2
b2
+
l2π2
c2
 )𝑡
 sen (
nπ
a
x) sen (
mπ
b
y) sen(
lπ
c
z) …(75) 
Aplicando las condiciones a la frontera: 
T(x, y, z, t) = ∑ ∑ ∑ Hnml
∞
l=1
∞
m=1
∞
𝑛=1 sen (
nπ
a
x) sen (
mπ
b
y) sen(
lπ
c
z) = ω(x, y, z) …(76) 
 
∫ ∫ ∫ sen (
nπ
a
x) sen (
mπ
b
y) sen (
lπ
c
z)
c
0
𝑤(𝑥, 𝑦, 𝑧)dxdydz 
𝑏
0
𝑎
0
= 
Hnml ∫ sen
2 (
nπ
a
x)
a
0
𝑑𝑥 ∫ sen2(
mπ
b
y)dy∫ sen2(
lπ
c
z)dz
c
0
b
0
 …(77) 
31 
 
Hnml = ∑ ∑ ∑
∫ ∫ ∫ sen(
nπ
a
x)sen(
mπ
b
y)sen(
lπ
c
z)
c
0
𝑤(𝑥,𝑦,𝑧)dxdydz 
𝑏
0
𝑎
0
 
∫ sen2(
nπ
a
x)
a
0
𝑑𝑥 ∫ sen2(
mπ
b
y)dy∫ sen2(
lπ
c
z)dz
c
0
b
0
∞
𝑙=1
∞
𝑚=1
∞
𝑛=1 …(78) 
 
T(x, y, z, t ) = ∑ ∑ ∑
∫ ∫ ∫ sen(
nπ
a
x)sen(
mπ
b
y)sen(
lπ
c
z)
c
0 𝑤(𝑥,𝑦,𝑧)dxdydz 
𝑏
0
𝑎
0 
∫ sen2(
nπ
a
x)
a
0 dx∫ sen
2(
mπ
b
y)dy
b
0 ∫ sen
2(
lπ
c
z)dz
c
0
∞
𝑙=1
∞
𝑚=1
∞
𝑛=1 𝑒
−α2(
n2π2
a2
+ 
n2π2
b2
+
l2π2
c2
)𝑡
 sen (
nπ
a
x) sen (
mπ
b
y) sen (
lπ
c
z) … (79) 
 
T(x, y, z, t ) = ∑ ∑ ∑
8
𝑎𝑏𝑐
∫ ∫ ∫ sen (
nπ
a
x) sen (
mπ
b
y) sen (
lπ
c
z)
c
0
𝑤(𝑥, 𝑦, 𝑧)dxdydz𝑒
−α2(
n2π2
a2
+ 
n2π2
b2
+
l2π2
c2
)𝑡
 sen (
nπ
a
x) sen(
mπ
b
y)sen(
lπ
c
z) 
𝑏
0
𝑎
0
∞
𝑙=1
∞
𝑚=1
∞
𝑛=1 …(80) 
 
Solución particular para la ecuación de transferencia de energía en una barra. 
 
 
 
 
 
 
 
32 
 
4.5 Caso V – Ecuación de difusión de masa 
 
Figura 5 - Placas 
Se obtendrá el desarrollo para obtener el perfil de concentraciones de la especie B, 
se propone una reacción química de primer orden entre la especie A y B para 
producir la especie C que será la capa de recubrimiento, a continuación se presenta 
dicha reacción química: 
 
𝐴 + 𝐵 
𝑘1
→ 𝐶 
La reacción propuesta es irreversible, equimolar y con una cinética de primer orden. 
Se consideran las siguientes suposiciones: 
 
I. Fluido newtoniano e incompresible. 
II. Proceso en estado estacionario. 
III. Las propiedades físicas de la fase fluida se mantienen constantes, incluyendo 
también al coeficiente difusión DAB. 
IV. Flujo laminar. 
Se retoman las suposiciones citadas anteriormente y se parte de la ecuación de 
continuidad para coordenadas rectangulares con densidad y coeficiente de difusión 
constante (Bird, 1976). 
 
 
De la ecuación, se observa que el efecto de la difusión, será mayor en el eje 𝑦 con 
respecto al eje 𝑥, debido a que en las abscisas, ésta dominada por el efecto 
conectivo, puesto que la velocidad del fluido se encuentra en esa dirección por ésta 
razón, en la ordenada, el efecto difusivo será mayor en comparación al eje x. 
 
33 
 
 
Por lo que la Ec. puede reducirse a: 
Vx
∂CB
∂x
 = DAB ∗
∂2CB
∂y2
 … (81) 
Partiendo de un perfil lineal, de acuerdo a la figura: 
Vx =
V
H
*y … (82) 
Sustituyendo el perfil en la Ec. Original, es decir en la Ec. 81: 
(
𝑉 
𝐻
𝑦)
∂CB
∂x
= DAB
∂2CB
∂y2
 … (83) 
𝜂 =
𝑦
𝐻
 ; 𝜑 = Vx
𝑉
 ; Ω = 𝐶𝐵
𝐶𝐵0
 ; 𝜉 = 𝑥
𝐿
 
Se adimensionaliza la Ec. 83 y la ecuación de difusión se reescribe como: 
𝜑 =
Vx
𝑉
= 𝜂 … (84) 
∂CB
∂x
=
∂CB
∂Ω
∂Ω
∂𝜉
 
∂𝜉
∂𝑥
 … (85) 
∂2CB
∂y2
 
∂
∂y∂𝜂
(
∂Ω
∂𝜂
) = 
𝐶𝐵0
H02
 
∂2Ω
∂𝜂2
 …(86) 
𝜂𝑈
𝐶𝐵0
𝐿
∂2CB
∂y2
= DAB
CB0
H02
∂2Ω
∂𝜂2
 … (87) 
∂Ω
∂𝜉
= 
DAB
𝑈
∗
L
H02
 ∂
2Ω
∂𝜂2
 … (88) 
Pe= U∗H0
2
L∗DAB
 … (88.a) 
 El Peclet de masa, como número adimensional, relaciona la velocidad de advección 
de un flujo y la velocidad de difusión. Es equivalente al número de Prandtl para difusión 
térmica. 
∂Ω
∂𝜉
=
1
𝑃𝑒
∂2Ω
∂𝜂2
 … (89) 
 Ω = 𝜓(𝜉)Γ(𝜂)=0 … (90) 
Se divide entre la ecuación 6 y se iguala a una constante de separación: 
𝜓′′(x)
𝜓(𝑥)
=
1
𝑃𝑒
1
𝜂
Γ′′(𝜂)
Γ(𝜂)
= −𝜆2 … (91) 
 𝜓’(𝜉)+λ2 𝜓 (𝜉)=0 … (92) 
34 
 
La solución es exponencial, tal como se ha observado anteriormente, para el caso de la 
primera ecuación diferencial parcial: 
Γ′′(𝜂) + 𝜆2𝑃𝑒𝜂 ∗ Γ = 0 … (93) 
Modelo de solución Airy con condiciones a la frontera: 
𝜓(𝜉) = 𝑐1𝑒
−𝜆2∗𝜉 … (94) 
𝑃(𝜂) = 𝐵𝐴𝑖𝑟𝑦𝐴𝑖 (
−𝑃𝑒∗𝜆2∗𝜂
(−𝑃𝑒∗𝜆2)
2
3
) + 𝐶𝐴𝑖𝑟𝑦𝐵𝑖 (
−𝑃𝑒∗𝜆2∗𝜂
(−𝑃𝑒∗𝜆2)
2
3
) …(95) 
La ecuación 95, muestra la solución particular de la ecuación de difusión de masa, para la 
distribución de concenctraciones: 
Ω(ξ, η) = 𝑐1𝑒
−𝜆2ξ [𝑐2𝐴𝑖𝑟𝑦𝐴𝑖 (
 −𝑃𝑒𝑀 𝜆
2η
(−𝑃𝑒𝑀 𝜆
2)
2
3
) + 𝑐3𝐴𝑖𝑟𝑦𝐵𝑖 (
−𝑃𝑒𝑀 𝜆
2η
(−𝑃𝑒𝑀 𝜆
2)
2
3
)] …(96) 
Ω(ξ, η) = 𝑒−𝜆
2ξ [𝑐4𝐴𝑖𝑟𝑦𝐴𝑖 (
 −𝑃𝑒𝑀 𝜆
2η
(−𝑃𝑒𝑀 𝜆
2)
2
3
) + 𝑐5𝐴𝑖𝑟𝑦𝐵𝑖 (
 −𝑃𝑒𝑀 𝜆
2η
(−𝑃𝑒𝑀 𝜆
2)
2
3
)] …(97) 
Las condiciones de frontera que se utilizarán para obtener las constantes de 
integración de la Ec. 97 son las siguientes: 
 Ω(𝜉, 0) = 𝑓𝑖𝑛𝑖𝑡𝑎 … (98) ; Ω(𝜉, 1) = 0 …(99) ; Ω(0, 𝜂) = 1 …(100) 
Al aplicar la condición a la frontera dada por la Ec. 98 a la Ec.97 se tiene que 𝑐5 = 0. 
Se aplica la condición de frontera dada por la Ec. 100. 
1 = ∑ 𝐶4
∞
𝑛=1
𝐴𝑖𝑟𝑦𝐴𝑖 (
 𝐷𝑎𝐼𝐼 − 𝑃𝑒𝑀 𝜆𝑛
2η
(−𝑃𝑒𝑀 𝜆𝑛
2)
2
3
)…(101) 
Para obtener la función peso de las funciones de Airy y poder aplicar el principio de 
ortogonalidad, se recurre a la estructura de la ecuación diferencial ordinaria, en este 
caso dada por la Ec. 93 y se obtiene que η es la función peso. 
1
2
∫ 𝐴𝑖𝑟𝑦𝐴𝑖 (
 − 𝑃𝑒𝑀 𝜆𝑛
2η
(−𝑃𝑒𝑀 𝜆𝑛2)
2
3
)𝑑η = 𝐶4
1
0
∫ [𝐴𝑖𝑟𝑦𝐴𝑖 (
−𝑃𝑒𝑀 𝜆𝑛
2η
(−𝑃𝑒𝑀 𝜆𝑛2)
2
3
)]
2
η 𝑑η
1
0
…(102) 
 
35 
 
Con lo que la constante quedan como: 
𝐶4 =
∫ 𝐴𝑖𝑟𝑦𝐴𝑖 (
 − 𝑃𝑒𝑀 𝜆𝑛
2η
(−𝑃𝑒𝑀 𝜆𝑛2)
2
3
)𝑑η
1
0
∫ [𝐴𝑖𝑟𝑦𝐴𝑖 (
 − 𝑃𝑒𝑀 𝜆𝑛2η
(−𝑃𝑒𝑀 𝜆𝑛2)
2
3
)]
2
η 𝑑η
1
0
 … (103) 
Finalmente se sustituye 𝑐4 en la siguiente ecuación y se tiene la distribución de 
concentraciones 
Ω(ξ, η) = ∑ 𝑒−𝜆𝑛
2ξ
∞
𝑛=1
[𝑐4𝐴𝑖𝑟𝑦𝐴𝑖 (
−𝑃𝑒𝑀 𝜆𝑛
2η
(−𝑃𝑒𝑀 𝜆𝑛
2)
2
3
)]… (104) 
4.6 Caso VI – Método perturbativo y capa límite 
 
Se resolverá la ecuación de Blasius, a través de una perturbación δ, método de expansión 
para mecánica cuántica. El método consiste en remplazar los términos no lineales, por 
una δ y después será tratada como un parámetro pequeño llamado capa límite. Las series 
de potencias para la delta convergen porque en las vecindades de la δ, no se presentan 
transiciones abruptas. (Kumar Datta, 2003) 
Ahora bien, se usará la perturbación δ, para la solución a la ecuación de Blasius, que 
describe nuevamente el perfil de velocidades, pero ésta vez sobre la capa límite, que se 
apreciar cuando un fluido corre sobre una placa o tubería. 
Para la ecuación: 
𝜕3𝑉
𝜕𝑥
+
𝜕2𝑉
𝜕𝑥
𝑉(𝑥) = 0… (105) 
Sujeto a las condiciones a la frontera 
V(0)=V´(0) …(106) ; V´(∞)=1 …(107) 
Nos interesa calcular el valor de Γ = V´´(0) del perfil de velocidades. Pese a que es 
complicado realizarse de manera computacional, tiene un rol importante en muchos 
análisis físicos. A través de métodos numéricos se tiene un valor muy aproximado a la 
cifra real con valor de Γ = V´´(0) = 0.46960. (Kumar Datta, 2003) 
 
 
36 
 
Se realizará el cálculo analítico a través del método de la delta perturbativa. Se hará a 
continuación cambiando una V por el parámetro delta, es decir: 
𝜕3𝑉
𝜕𝑥
+
𝜕2𝑉
𝜕𝑥
(𝑉(𝑥))
𝛿
= 0… (108) 
Si se asignará un valor de uno a la δ, se tendría la forma original de Blasius. Para resolver 
por el método perturbativo, se expresará la velocidad V(x) en series de potencias de δ es 
decir: 
𝑉 = 𝑉0(𝑥) + 𝛿𝑉1(𝑥) + 𝛿
2𝑉2(𝑥) + 𝛿
3𝑉3(𝑥)…+ 𝛿
𝑛𝑉𝑛(𝑥)… (109) 
Sustituyendo la expansión en la ecuación tenemos una secuencia de ecuaciones lineales, 
asociadas a condiciones a la frontera: 
𝜕3𝑉0
𝜕𝑥
+
𝜕2𝑉0
𝜕𝑥
= 0… (110) 
𝑉0(0) = 𝑉0´(0) = 0, 𝑉1´(∞) = 1 
𝜕3𝑉1
𝜕𝑥
+
𝜕2𝑉1
𝜕𝑥
= −
𝜕2𝑉0
𝜕𝑥
ln(𝑉0) = 0… (111) 
𝑉1(0) = 𝑉1´(0) = 𝑉1´(∞) = 0 
𝜕3𝑉2
𝜕𝑥
+
𝜕2𝑉2
𝜕𝑥
= −
𝜕2𝑉1
𝜕𝑥
ln(𝑉0) −
𝜕2𝑉0
𝜕𝑥
(
𝑉1
𝑉0
+
1
2
𝑙𝑛2(𝑉0))… (112) 
𝑉2(0) = 𝑉2´(0) = 𝑉2´(∞) = 0 
Algo típico en las series de expansión por método perturbativo, es que la primera 
ecuación generada es siempre homogénea, mientras que todas las demás restantes no lo 
son teniendo todas sus límites o condiciones a la frontera definidos. 
Ahora bien integrando la ecuación diferencial (110) de orden cero, arroja una solución 
exponencial, que puede ser calculada como: 
𝑑3(𝑒𝜆𝑥)
𝑑𝑥3
+
𝑑2(𝑒𝜆𝑥)
𝑑𝑥2
= 0… (113) 
Sustituyendo: 𝑑
2(𝑒𝜆𝑥)
𝑑𝑥2
= 𝜆2𝑥 y lo mismo por la parte cúbica es decir 𝜆3𝑥 
Resta obtener los valores propios como: 
(𝜆3 + 𝜆2) 𝑒𝜆𝑥 = 0… (114) 
Al obtener las raíces factorizando, se sabe que: 𝜆1 = −1, 𝜆2,3 = 0. 
 
37 
 
Solución general para el orden cero queda descrito como: 
𝑉0(𝑥) = 𝑒
−𝑥 + 𝑥 − 1… (115) 
Para la condición, cuando esté valuado en 0, da como solución: 
𝛤 = 𝑉0´´(0) = 1… (116) 
Valor que ésta bastante próximo, a la solución real, sin embargo, si se evalúa el segundo 
orden, el acercamiento es aún mayor. 
Solución de primer orden 
La segunda ecuación presenta la ventaja de que para la segunda derivada de 𝑉1(𝑥) para 
x = 0, a eso le sigue que 𝑉1(𝑥)´´ pueda ser obtenido a través de la integral: 
𝑉1(0)´´ = ∫ 𝑒
−𝑦𝑙𝑜𝑔(𝑒−𝑦 + 𝑦 − 1)
∞
0
dy… (117) 
Que no presenta una solución analítica conocida, sin embargo a través de una integración 
numérica se obtiene que: 
𝑉1(0)´´ = −2.1333275 … (118) 
Que está obviamente muy lejos del resultado esperado. Usando una aproximación de 
Pade, que supone la solución de la derivada inmediata anterior, donde se tiene que variar 
la delta, de tal manera que el valor se encuentre entre 0 y 1, para que converja a 0.46960, 
más la solución del recién calculado por el parámetro delta, nos acerca mucho más a la 
solución exacta: 
𝛤 = 1 − 2.333275 ∗ 𝛿 = 0.31914 … (119) 
Se puede apreciar que el valor es por mucho más cercano al dela primera derivada, sin 
embargo, a mayor extensión de la serie el error esperado será aún menor, es decir a 
mayor orden de perturbación, mayor será la exactitud de nuestra solución. 
Solución de segundo orden 
Los cálculos para ésta ecuación se obtienen de forma directa, sin embargo requieren de 
un poco más de álgebra que para el cálculo de orden 0, puesto que se requiere de cambio 
de variable para la resolución de la integral. La solución a la ecuación de segundo orden 
ésta dada por: 
𝑉2´´(0) = ∫ 𝑒
−𝑦
∞
0
𝑦𝑙𝑜𝑔(𝑉0)𝑑𝑦 +
1
2
∫ 𝑒−𝑦
∞
0
(1 − 𝑦)𝑙𝑜𝑔2𝑉0(𝑦)𝑑𝑦… (120) 
Por cambio de variable, se llega a 
𝑉2´´(0) = ∫
𝑒−𝑦
𝑉0(𝑦)
∞
0
𝑑𝑦 +
1
2
∫ 𝑒−𝑧
𝑦
0
(𝑧)𝑙𝑜𝑔𝑉0(𝑧)𝑑𝑧… (121) 
38 
 
Valuado en 0 
𝑉2´´(0) = 
1
2
((𝑉1(0)´´))
2
− 𝑉1´(0) +
1
2
= 5.83101 … (122) 
Se usa nuevamente la técnica de Pade para aproximar: 
𝛤 = 1 − 2.333275 ∗ 𝛿 + 5.83101 ∗ 𝛿2 = 0.42901 … (123) 
La solución obtenida para el segundo orden tiene un error del 14%, lo cual ya se puede 
considerar correcto, sin embargo, si se quisiera desarrollar una última serie para ver el 
nuevo cambio, el camino a seguir debe ser el mismo. Para el cálculo iterativo, se requiere 
de la solución del segundo orden, más la integración del de tercer orden, que como puede 
suponerse nos arrojará un término muy pequeño, ya que desde el segundo orden el error 
presenta un porcentaje pequeño. 
5 Método del elemento finito (MEF ó FEM) 
 
La técnica del elemento finito, es hoy en día el segundo método numérico más usado en 
el mundo, después de la transformada de Fourier, no sólo por Ingenieros sino también 
para los físicos, por su gran aplicación en la mecánica y electromagnetismo. Como se 
estudió en los antecedentes, éste método tiene sus raíces en los trabajos de residuos 
ponderados de Petrov Galerkin y en los métodos variacionales de R. Ritz en los espacios 
de Sobolev. 
Explicar el método por su extensión y complejidad, requiere de un curso completo y 
debido a que el tema es para este proyecto, sólo un medio para la solución de ecuaciones 
diferenciales, no se desarrollará tan a fondo. 
Ahora bien el punto de inicio para dicho método, comienza a partir de la siguiente técnica, 
conocida como la solución débil del sistema (Débil porque al final del método se tiene un 
error muy cercano al cero) sin embargo dicho error conlleva a que el método haya perdido 
exactitud, de acuerdo a lo que se explica en el siguiente apartado. 
5.1 Solución débil 
 
Sea V en un espacio real R, y se tiene una función lineal en V es decir: 
 l.v R Э 
1.-) l(cV) = cl(v) ∀ Cϵ R Λ ∀v ϵ V 
2.-) l(V) = l(u) ∀ u ϵ V 
 
39 
 
Para una función del tipo: 
(𝑓 + 𝑔)(𝑥) = 𝑓(𝑥) + 𝑔(𝑥) ∀ f ϵ R 
Se requiere llegar a la ecuación general de Poisson, que es derivada de la suposición 
anterior, aplicada a un sistema dinámico, podría verse como el siguiente ejemplo: 
 
Figura 6 - Banda elástica 
Donde se puede apreciar que se trata de una fuerza aplicadahacia una dirección (En éste 
particular caso a la derecha) a una presión p(x) donde la x es la variable fuerza, en una 
dirección con longitud de flecha l. Este problema genera una ecuación de la siguiente 
forma: 
𝑑
𝑑𝑥
(𝐸𝐴
𝑑𝑢
𝑑𝑥
) + 𝑐𝑢 = 𝑝… (124) 
Con condiciones a la frontera BC1: u(l)-u=0, BC2: −EA du(0)
dx
+ F0 = 0 
Ahora bien, dado que se trata de una ecuación diferencial, sería más facil aproximarla por 
un método numérico, si se multiplicase por una “Test function’’ o en español una funcional 
‘’v’’ que se encargara de multiplicar nuestra función por una funcional tal, que nos permita 
más adelante discretizar nuestro sistema. (Tecnische Universität München, 2015) 
Si para el problema anterior se considera 𝑣(𝑙) = 0 y se multiplica la ecuación diferencial 
por la condición a la frontera 2, entonces nos queda: 
(−
𝑑
𝑑𝑥
(𝐸𝐴
𝑑𝑢
𝑑𝑥
) + 𝑐𝑢 − 𝑝 = 0 ) ∗ 𝑣 … (125) 
Usando la condición de la frontera 1, es decir cuando u(l) vale uno, recordando que u es 
la función desplazamiento, la ecuación 125 toma la forma: 
−𝐸𝐴
𝑑𝑢(0)
𝑑𝑥
+ 𝐹0 = 0 …(126) 
E integrando sobre la barra, que recorre una línea desde cero a la longitud de la misma 
(l): 
∫ (−
𝑑
𝑑𝑥
(𝐸𝐴
𝑑𝑢
𝑑𝑥
) + 𝑐𝑢 − 𝑝) ∗ 𝑣 𝑑𝑥 = 0… (127)
𝑙
0
 (Tecnische Universität München, 2015) 
 
40 
 
Integración por partes 
∫ 𝐸𝐴
𝑑𝑢
𝑑𝑥
∗
𝑑𝑣
𝑑𝑥
+ 𝑐𝑢𝑣 − 𝑝𝑣 𝑑𝑥 … (128)
𝑙
0
 
Es éste el punto principal del método, pues en ello recae la exactitud de la solución al 
método numérico propuesto (MEF) es decir que se tiene que encontrar una función 
espacial U o de desplazamiento, tal que para cualquier función v, la ecuación quede 
descrita de la siguiente forma: 
∫ 𝐸𝐴
𝑑𝑢
𝑑𝑥
∗
𝑑𝑣
𝑑𝑥
+ 𝑐𝑢𝑣 𝑑𝑥 = ∫ 𝑝𝑣 𝑑𝑥 − 𝐹0𝑣(0)… (129)
𝑙
0
𝑙
0
 
B(u,v) = F(v) 
Observar que se adquiere una forma bilineal. 
Bajo condiciones suficientes, la solución de la ecuación diferencial fuerte, será idéntica a 
la solución diferencial débil. 
Se sabe que es bilineal pues para toda: 
𝐵(𝑢, 𝛼𝑣1 + 𝛽𝑣2) = 𝛼𝐵(𝑢, 𝑣1) + 𝛽𝐵(𝑢, 𝑣2)… (130) 
𝐵(𝛼𝑢1 + 𝛽𝑢2, 𝑣) = 𝛼𝐵(𝑢1, 𝑣) + 𝛽𝐵(𝑢2, 𝑣)… (131) 
Que además cumple condiciones de simetría, es decir: 
𝐵(𝑢, 𝑣) = 𝐵(𝑣, 𝑢)para toda u, v y 𝐵(𝑢, 𝑣) > 0 es positiva definida para u≠ 0 
Es importante recordar de aquí en adelante, que se usaron los trabajos de Ritz sobre 
principios variacionales, puesto que aquí es puramente algebra ya que entraremos en 
espacios vectoriales de Rn, por ello mismo se enfatizarán algunas propiedades para las 
funciones espacio, en espacios vectoriales: 
- La suma de dos elementos de un conjunto, es un elemento del conjunto. 
- El producto de un escalar por un elemento del conjunto, es y será un elemento del 
conjunto. 
Recordar, que para funciones continuas: 
A un subespacio de (E) se dice que produce (spann) un conjunto de funciones en el 
espacio V, si cada función v ϵ V, entonces puede expresarse: 
𝑣 = ∑𝑎𝑖𝐸𝑖 donde ai es un escalar … (132) 
Donde se puede apreciar que la sumatoria tendría un infinito número de elementos, es 
decir, que es infinitamente dimensional la matriz. También se observa que E, soporta 
cualquier función v, siempre y cuando la B y F0 sean linealmente independientes. 
41 
 
5.2 Acercamiento de Galerkin 
 
Establece que no se deben de probar todas las posibles funcionales, sino sólo aquellas 
(número finito) que sean linealmente independientes a través de las funciones de forma 
(Ver antecedentes) es decir para la 𝑁𝑗 , 𝑑𝑜𝑛𝑑𝑒 𝑗 = 1,2, … ,𝑚. 
En segundo término, se dice que se aproxime la función U o de desplazamiento, como la 
sumatoria de los términos linealmente independientes es decir: 
𝑈 = ∑ 𝑎𝑖𝑁𝑖
𝑚
𝑖=1 …(133) Donde 𝑎𝑖 son coeficientes. (Krasnov, 1976) 
Entonces se puede resumir la forma débil de la ecuación diferencial parcial de cualquier 
orden, podrá ser escrita cuando se encuentren los coeficientes 𝑎𝑖 , i=1, 2,…, m de 𝑈 =
∑ 𝑎𝑖𝑁𝑖
𝑚
𝑖=1 tal que se cumpla a condición a la frontera 1, es decir que: 
𝐵(𝑈,𝑁𝑗) = 𝐵 (∑𝑎𝑖𝑁𝑖,𝑁𝑗
𝑚
𝑖=1
) = 𝐹(𝑁𝑗)… (134) 
∑𝑎𝑖𝐵(𝑁𝑖,𝑁𝑗)
𝑚
𝑖=1
= 𝐹(𝑁𝑗) 𝑝𝑎𝑟𝑎 𝑡𝑜𝑑𝑎 𝑗 = 1,2, … ,𝑚. 
De manera general se puede decir que el plan de trabajo para elemento finito es el 
siguiente (Tomado de un resumen de elemento finito computacional de la TUM) 
 
Figura 7 - Metodología FEM (Tecnische Universität München, 2015) 
42 
 
Este esquema representa de manera muy sencilla, la metodología a seguir paso a paso, 
para el cálculo del (FEM- (Nevena, 2015)) ya sea de manera analítica o aproximada a 
través de un software de computadora (En éste caso será Matlab) Mismo que será 
desarrollado más adelante, después de abordar el cálculo del problema 1D, es decir el 
más sencillo que se ha estado llevando a cabo, el problema de la barra. 
 
Figura 8- Barra Segmentada (Tecnische Universität München, 2015) 
Supongamos que se analiza tal problema, estudiado anteriormente, para la obtención de 
la ecuación de Poisson, se puede apreciar que la misma ha sido seccionada en una 
infinidad de tramos, o elementos, es decir que se tienen e-n elementos en donde la 0 
=x1<x2<x3<…xm= l, en un dominio 𝛺𝑒 = 𝑥𝑘 < 𝑥 < 𝑥𝑘+1. 
Dado que la barra total, se ha repartido en muchos partes que pueden ser o no de 
diferente tamaño, se parte de la suposición de que los elementos son del mismo tamaño y 
entonces pueden considerarse coordenadas locales (ζ) 
Donde x puede ser expresado de la siguiente manera: 
𝑥 = 𝑄𝑒 = 
1−𝜁
2
𝑥𝑘 + 
1+𝜁
2
𝑥𝑘 + 1…(135) (Nevena, 2015) 
Misma expresión, que se puede analizar y obtener de segmentar cada uno de los 
elementos, dentro de nuestra barra y obteniendo las funciones de forma N, para el caso 
1D, como se muestra a continuación (Líneas rojas) 
 
 
Figura 9- División de secciones 
 
 
43 
 
Donde Qe es la transformación de coordenadas globales en locales, a través de las 
funciones N, es decir: 
 
Figura 10- Funciones de forma 
 N1 = 1−𝜁
2
𝑥𝑘 y N2 = 1+𝜁
2
𝑥𝑘+1 
Éstas funciones son de gran importancia, pues la aproximación final (U) es obtenida a 
través de la combinación lineal de las mismas, es decir cuando se incorporan aquellas 
soluciones en la sumatoria, que son no lineales, como se muestra en la figura 11. 
 
Figura 11- Solución desplazamiento 
Una vez seleccionado el mallado que queremos realizar (En éste caso se trata de una 
barra o un elemento 1D, el cual no es mallado) ya se considera forma débil o discretizada 
de nuestra ecuación diferencial, posteriormente se procede a realizar un sistema de 
ecuaciones algebraicas a través de una matriz de elementos, que será construida por los 
métodos descritos anteriormente. No olvidar que para cualquiera que sea el orden de la 
EDP, el método será el mismo, sin embargo, para las de segundo orden dependiendo la 
geometría, las matrices de fuerza y la de elementos cambiarán significativamente si se 
trata de un elemento triangular o cuadrático. 
44 
 
5.3 Matriz de elementos (Element stiffness matrix) para el caso 1D 
 
Es el más sencillo de los problemas, ya que la matriz presenta una forma diagonal, puesto 
que las coordenadas son consecutivas, es decir, el elemento uno se junta con el dos, en 
un par de coordenadas, como también el dos con el tres y así sucesivamente. 
La matriz K o bien matriz de elementos, está definida como: 
𝑘𝑖𝑗 = 𝐵(𝑁𝑖, 𝑁𝑗) que es lo mismo a 𝐵(𝑁𝑖, 𝑁𝑗)= ∫ 𝐸𝐴
𝑑𝑁𝑖
𝑑𝑥
𝑑𝑁𝑗
𝑑𝑥
𝑙
0
𝑑𝑥 … (136) 
Que es igual a escribir: 
∑ ∫ 𝐸𝐴
𝑑𝑁𝑖𝑒
𝑑𝑥
𝑑𝑁𝑗𝑒
𝑑𝑥
𝑙
0
𝑛−1
𝑒=1
𝑑𝑥 …(137) 
Darse cuenta que si el resultado de la sumatoria es igual a cero, o bien Ni ó Nj son cero. 
Considérese ahora, que ninguna de las funciones de forma es igual a cero y que además 
los elementos del dominio 𝛺𝑒 con nodos en i,j, definen la matriz elemento para el caso 1D 
como: 
Ke = (
𝑘𝑖𝑖
𝑒 𝑘𝑖𝑗
𝑒
𝑘𝑗𝑖
𝑒 𝑘𝑗𝑗
𝑒 )… (138) 
Donde kij, es obtenida a través de la suma de todas