Ecuación del calor, equipo LUA
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación del calor. Grupo LUA |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Luis Carreras Hoyos, Lucía Gil Ruiz y Alejandra Hernández Sieber |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Conceptos previos
- 3 Modelización de la ecuación del calor con una dimensión
- 4 Cambio de coeficientes
- 5 Cambio de condiciones
- 6 Cambios de condiciones de frontera
- 7 Solución fundamental de la ecuación del calor
- 8 Referencias
1 Introducción
La ecuación del calor es una de las ecuaciones diferenciales parciales más fundamentales en la física y la ingeniería, utilizada para modelar la difusión del calor en sólidos, líquidos y gases. Su estudio proporciona una comprensión profunda de los procesos de transferencia de calor y la evolución temporal de la temperatura en sistemas físicos. En este trabajo, nos enfocamos inicialmente en la ecuación del calor en una dimensión, explorando su formulación matemática y calculando soluciones para diversas condiciones iniciales y de frontera, además de examinar las implicaciones físicas de estas soluciones. En concreto, estudiaremos la variación temporal de la temperatura en función de la posición a lo largo de una varilla metálica, considerando condiciones iniciales y de frontera diversas, como condiciones de contorno de tipo Dirichlet y Neumann. Posteriormente, ampliaremos nuestro enfoque para abordar la ecuación del calor en dos dimensiones. Aunque este estudio será menos extenso en comparación con el caso unidimensional, nos permitirá explorar cómo se extienden los conceptos y métodos desarrollados previamente a sistemas bidimensionales.
2 Conceptos previos
Antes de comenzar con el trabajo, en esta sección vamos a definir algunos conceptos previos necesarios para la completa comprensión del mismo.
2.1 Contextualización y deducción de la ecuación del calor
Procedemos a deducir la ecuación del calor. En primer lugar, la Ley de enfriamiento de Newton lleva el nombre de Isaac Newton, quien la formuló en el siglo XVII. Esta es una ley empírica cuya idea es que la velocidad a la que disminuye (o aumenta) la temperatura [math] T(t) [/math] de un cuerpo es proporcional a la diferencia entre la temperatura ambiental [math] T_a [/math] y la del cuerpo. Es decir,
Posteriormente, J.Fourier estableció que dichas temperaturas no eran constantes. Consideró [math] V [/math] un volumen de control y definió [math] e(\mathbf{x},t) [/math] como la densidad de energía calorífica, siendo esta la cantidad de energía térmica por unidad de volumen en un medio, con unidades [math] \frac{cal}{m^3} [/math].
Dado que el calor es una forma de energía, es natural hacer uso de la ley de conservación de la energía durante la deducción de la ecuación del calor. Esta postula que sobre [math] V [/math], la variación de la energía calorífica se debe al balance entre el calor que entra y el que sale, más la consideración de una posible producción externa. Por tanto, se deduce:
Siendo [math] \overrightarrow{q} [/math] el flujo de calor, medido en [math] \frac{cal}{m^2s} [/math] y [math] r [/math] una posible fuente externa de calor, medida en [math] \frac{cal}{m^3s} [/math].
A continuación, por el Teorema de la divergencia se tiene que sea [math] \Omega \in \mathbb{R^3} [/math] un dominio Lipschitz, entonces:
Aplicando a la Ley de conservación de la energía el Teorema de la divergencia, obtenemos
Como se cumple para todo volumen [math] V \in \Omega[/math], el integrando debe ser nulo, llegando a que [math] \frac{\partial e(\mathbf{x},t)}{\partial t} + div(\overrightarrow{q})=r [/math].
Para continuar, asumimos lo siguiente:
- Ley de Fourier de la conducción del calor: el flujo de calor [math] \overrightarrow{q} [/math] a través de un material es directamente proporcional al producto del gradiente de temperatura [math]T(\mathbf{x},t)[/math] y la conductividad térmica [math] k [/math] del material. En términos de fórmulas, sería [math] \overrightarrow{q} = -k \cdot \nabla T(\mathbf{x},t) [/math]
- La energía térmica es proporcional a la temperatura, esto es [math] e(\mathbf{x},t) = c \cdot T(\mathbf{x},t) [/math], siendo [math] c [/math] la constante de calor específico.
Con ello, hemos llegado a que
Donde en la última igualdad hemos empleado que [math] div(\nabla F) = \Delta F [/math]. Agrupando términos, se concluye que la ecuación del calor </math> sería:
2.2 Constante de conductividad térmica y calor específico
A continuación, introducimos tres definiciones:
- Constante de conductividad térmica: La constante de conductividad térmica, a la que hemos denotado como [math] k [/math], es una propiedad intrínseca de los materiales la cual describe su capacidad para conducir el calor y se mide en [math] \frac{cal}{ºC s m}[/math]. Cuanto mayor sea la constante de conductividad térmica, mayor será la capacidad del material para conducir el calor.
- Calor específico: El calor específico, denotado previamente como [math] c [/math], es una propiedad de cada material que indica cuánto calor es necesario para elevar la temperatura de una unidad de masa de esa sustancia en una unidad de temperatura, medida [math] \frac{cal}{ºC g}[/math]. Cuanto mayor sea el calor específico de una sustancia, más calor se requiere para aumentar su temperatura.
- Constante de difusión térmica: La constante de difusión térmica, a la que hemos denotado como [math] D [/math], es el cociente de la constante de conductividad térmica entre el calor específico.
2.3 Principio del máximo
El principio del máximo para la ecuación del calor establece que si la función [math]u(x,t) \in C^{2,1}(Q_{T}) \cap C(\overline{Q}_{T}) [/math] cumple la condición [math]u_{t}-D\Delta u \leq 0 [/math] para [math](x,t) \in Q_{T}[/math]. Entonces [math] u(x,t) [/math] verifica [math] \max_{x\in \partial_p Q_{T}} u(x,t)= \max_{x\in \overline{Q_{T}}} u(x,t) [/math].
Donde [math]Q_{T}= \Omega \times (0,T) [/math], siendo [math] \Omega \subset R^{3}[/math].
2.4 Convolución
La convolución es una operación matemática que combina dos funciones para describir la superposición entre ambas. La convolución de las funciones [math] f [/math] y [math] g [/math] se denota por [math] f *g [/math] y se define por:
En este trabajo, usaremos esta operación para obtener una solución de la ecuación del calor dado un dato inicial concreto, como veremos posteriormente.
3 Modelización de la ecuación del calor con una dimensión
Para estudiar el comportamiento de la ecuación del calor, vamos a comenzar tratando el siguiente problema:
Se considera una varilla metálica que ocupa el intervalo [math] [0,1] [/math] y que se encuentra aislada por su superficie lateral, de manera que la conducción de calor sólo se produce en la dirección longitudinal. La temperatura inicial de la varilla es [math] 0 °C [/math]. En el extremo izquierdo se consigue mantener la temperatura a [math] 0 °C[/math] mientras que en el derecho la temperatura es siempre de [math] 1°C [/math].
Aplicando la fórmula de la ecuación del calor, definiendo [math] u(x,t) [/math] como la función que determina la temperatura y tomando como constante de conductividad térmica [math] k=1~\frac{cal}{ºC s m} [/math] y calor específico [math] c=1~\frac{cal}{ºC g} [/math], modelizamos el comportamiento de la temperatura mediante el correspondiente sistema de EDP: [math] \begin{cases} u_t-u_{xx}=0 & x\in(0,1),t\gt0\\ u(x,0)=0&x\in(0,1)\\ u(0,t)=0&t\gt0\\ u(1,t)=1&t\gt0 \end{cases}[/math]
En primer lugar, cuando [math] t \rightarrow \infty[/math], [math] u_t \sim 0 [/math] y por ello, la solución no varía en tiempo, es decir, [math] u(t,x) \sim v(x) [/math], siendo [math] v(x) [/math] la solución estacionaria y cumpliendo que:
[math] \begin{cases} v_{xx}=0 & x\in(0,1)\\ v(0)=0\\ v(1)=1\\ \end{cases}[/math]
Resolviendo el sistema obtenemos [math] v(x)=x [/math], cuya representación es:
clear all
close all
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) x;
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
numero_puntos=100;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
X=linspace(extremo_izquierdo,extremo_derecho,numero_puntos); % vector del mallado espacial
plot(X,sol_estacionaria(X))
xlabel('Espacio');
ylabel('Temperatura');
title("Solución estacionaria")
Esta solución estacionaria se resta a la solución original obteniendo el cambio de variable [math] w(x,t)=u(x,t)-v(x) [/math]. Mediante este cambio, obtenemos nuestro sistema de EDP’s homogeneizado:
[math] \begin{cases} w_t-w_{xx}=0 & x\in(0,1),t\gt0\\ w(x,0)=-x&x\in(0,1)\\ w(0,t)=0&t\gt0\\ w(1,t)=0&t\gt0 \end{cases}[/math]
Tomamos ahora el cambio [math] z(x,t)=-w(x,t) [/math], donde obtenemos el sistema final: [math] \begin{cases} z_t-z_{xx}=0 & x\in(0,1),t\gt0\\ z(x,0)=x&x\in(0,1)\\ z(0,t)=0&t\gt0\\ z(1,t)=0&t\gt0 \end{cases} [/math]
Resolviendo el sistema por separación de variables, obtenemos como solución [math] z(x,t)=\sum_{k=1}^\infty\frac{2(-1)^{k+1}}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t} [/math]. Deshaciendo los cambios que hemos realizado, obtenemos la solución al problema original, [math] u(x,t)=\sum_{k=1}^\infty(-1)\frac{2(-1)^{k+1}}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t}+x [/math]. Representamos la solución [math] u(x,t) [/math] para [math]x \in [0,1] [/math] en el intervalo de tiempo [math]t \in [0,1] [/math] tomando solo los 10 primeros términos de la serie.
En la imagen, podemos ver como a medida que agregamos términos, la solución se vuelve más precisa. Además, vamos a comprobar como las condiciones iniciales de nuestro sistema concuerdan con la solución obtenida. En la posición [math] x=0 [/math] la solución es la constante [math] 0 [/math] y en la posición [math] x=1[/math] la solución es constantemente [math] 1 [/math] . Sin embargo, a medida que aumenta la posición [math] x [/math] en el tiempo inicial, la solución no vale constantemente 0, sino que oscila temporalmente. Esto es debido a que la solución original no es una función continua cerca de [math](x,t)=(1,0) [/math]. Por un lado, tenemos la condición [math]u(x,0)=0[/math], mientras que por el otro lado [math]u(1,t)=1[/math]. Debido a la incapacidad de las series de Fourier para representar de manera precisa las discontinuidades, las sumas parciales de Fourier tienden a "sobrepasar" la discontinuidad en un intento de ajustarse a la función. Este suceso es conocido como el fenómeno de Gibbs y es el responsable de las fluctuaciones mencionadas.
A continuación, aumentamos el valor de la serie hasta [math]k=1000[/math] para ver cómo varia:
Como se puede apreciar en esta imagen, cuando [math] k=1000 [/math], este último caso no se aprecian estas oscilaciones. A medida que se agregan más términos a la serie, las oscilaciones disminuyen en amplitud.
Además, esta discontinuidad en el instante inicial provoca que la función no sea diferenciable en [math] t=0 ~s [/math]. De hecho, si calculamos la derivada espacial en la solución término a término de manera analítica en dicho instante, obtenemos que [math] \sum_{k=1}^{\infty} \frac{d}{dx} u_{k}(x,0) + \frac{d}{dx} x=1+\sum_{k=1}^{\infty} 2(-1)^{k} cos(k\pi x) [/math]. En particular, si evaluamos en [math] x=0 [/math]
Esta es una serie que no converge, pues toma los valores [math] -1 [/math] y [math] 1 [/math] alternadamente. De manera equivalente, evaluando en [math] x=1 [/math]
Esta última serie tampoco converge, pues el límite de sus términos es claramente no nulo.
Sin embargo, para instantes de tiempo no nulos la solución converge uniformemente y, por ello, se verifica que
Por lo tanto, concluimos que no tiene sentido analítico calcular la derivada espacial en [math] t=0 [/math].
Asimismo, se ve con claridad como la solución obtenida se aproxima al estado estacionario, pues converge a la función [math] v(x)=x [/math] cuando el tiempo es suficientemente grande.
3.1 Flujo del calor
Ahora vamos a analizar cómo varía el flujo del calor entrante y saliente en los extremos de la varilla a lo largo del tiempo. Para ello, conocemos que el flujo sigue la ecuación [math] \overrightarrow{q} = - \nabla u(x,t) [/math], teniendo en cuenta que el coeficiente de difusión vale [math] 1 [/math] en este caso. Considerando la solución con [math] 10 [/math] y [math] 1000 [/math] términos y evaluando el flujo en Matlab para ambos extremos, obtenemos las siguientes gráficas:
Recordamos que en el extremo izquierdo, la temperatura vale constantemente [math] 0°C [/math], mientras que en el extremo derecho queda fijada a [math] 1°C [/math]. Por ello, a medida que va transcurriendo el tiempo, el calor se difunde desde donde más temperatura hay a donde menos, es decir, del extremo izquierdo de la varilla al derecho.
Como en el instante inicial, toda la varilla se encuentra a temperatura [math] 0°C [/math], excepto el extremo derecho, el flujo del extremo izquierdo para instantes próximos al inicial debe valer [math] 0 ~ \frac{cal}{m^2 s} [/math] aproximadamente pues los puntos próximos a él tienen la misma temperatura. Tal y como se cumple cuando aumentamos la precisión para k=1000. A medida que transcurre el tiempo, el valor absoluto de este flujo irá aumentando pues los puntos próximos al extremo izquierdo irán adquiriendo temperatura. Sin embargo, el flujo en el extremo derecho tiene un comportamiento totalmente contrario. En el instante inicial, hay una discontinuidad abrupta de temperatura en el extremo derecho, lo que provoca un flujo cuyo valor absoluto es realmente elevado. Además, el calor se irá transmitiendo lentamente a lo largo de la varilla, provocando que cada vez el valor absoluto del flujo sea más pequeño.
Por otra parte, cabe destacar que, ambos flujos toman valores negativos, lo cual indica que su sentido es opuesto al propio eje direccional, provocando que sea saliente en el extremo izquierdo y entrante en el derecho.
Cabe destacar que, tal y como hemos mencionado antes, la solución no es derivable con respecto al espacio en el instante inicial, y por ello, no tiene sentido estudiar el flujo en [math] t=0 [/math].
Finalmente, es imprescindible mencionar que esta gráfica es excelente a la hora de estudiar la convergencia de la solución a la solución estacionaria. El instante de tiempo en el que ambos flujos se estabilizan tomando el mismo valor simboliza que, en la varilla, el calor que entra por un extremo es el mismo que el que sale por el otro extremo, lo que obviamente es una consecuencia de haber alcanzado el estado estacionario.
3.2 Código de Matlab
Todas las gráficas que se han representado en esta sección, son resultado del siguiente código:
clear all
close all
% COEFICIENTES DE FOURIER
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t));
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) x;
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.01;
punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=1;
paso_tiempo=0.0001;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
w=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
k_end=1000; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
for ks=1:k_end % queremos resolver el problema para distintos valores de k
w=w+w_k(X,T,ks,coeficiente_conductividad_termica/calor_especifico,c_k(X(1,:),ks));
u=w+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
if ks <= 10 % representamos las gráficas de a lo sumo los 10 primeros términos de la solución
figure(1)
sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
subplot(2,min(ceil(k_end/2),5),ks)
hold on
mesh(X,T,u)
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
colormap('jet')
colorbar
end
% COMPARACIÓN DE LA SOLUCIÓN PARA k=10 Y k=k_end
if (k_end>10) & (ks==10)
figure(2)
sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
subplot(1,2,1)
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
title("k = "+num2str(ks))
figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
sgtitle('Representación del flujo con respecto al tiempo en los extremos para distintos valores de k')
subplot(1,2,1)
deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
plot(T(:,1),flujo0,'r','linewidth',1.5)
hold on
plot(T(:,1),flujo1,'b','linewidth',1.5)
xlabel('Tiempo')
ylabel('Flujo')
legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
title("k="+num2str(ks))
elseif (ks==k_end)
figure(2)
if (ks>10)
subplot(1,2,2)
end
mesh(X,T,u)
colorbar
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
if (ks>10)
subplot(1,2,2)
end
deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
plot(T(:,1),flujo0,'r','linewidth',1.5)
hold on
plot(T(:,1),flujo1,'b','linewidth',1.5)
axis([0,1,-100,0])
xlabel('Tiempo')
ylabel('Flujo')
legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
title("k="+num2str(ks))
end
end
En el código, lo primero que hemos hecho ha sido definir los coeficientes de Fourier, c_k, y el término k-ésimo de la solución del problema después de haber restado la solución estacionaria,w_k. Además, se define la solución estacionaria del problema, sol_estacionaria.
Luego hemos establecido valores para los coeficientes de conductividad térmica y el calor específico y se ha discretizado tanto el tiempo como el espacio para crear una malla espaciotemporal X y T mediante la función meshgrid. De esta manera, pasamos de un problema continuo a un problema discreto con un conjunto finito de puntos. Esta discretización nos permite representar la función gráficamente y utilizar métodos numéricos para el cálculo de los coeficientes de Fourier, como veremos posteriormente.
A continuación, calculamos la solución w homogeneizada con k-ésimos términos [math] (k=1,\dots,100) [/math], mostrando en gráficas las distintas soluciones y su evolución con respecto al espacio y tiempo. Además, graficamos la solución para [math] k=10 [/math] y [math] k=1000 [/math] y calculamos los flujos en los extremos, empleando la función gradient. Aunque analíticamente no es posible calcular dichos flujos en el instante inicial, hemos hecho uso de esta función para calcular la derivada de manera numérica. Esto permite que el lector no se vea obligado a calcular dicha derivada por sí mismo, pero puede cometer errores como esta aproximación, así siempre debemos corroborar los resultados obtenidos en estas gráficas con la base teórica.
Cabe destacar que todos estos datos están al comienzo del programa con la idea de que sean fácilmente identificables y modificables para el lector, y así poder estudiar distintas situaciones iniciales que no abordamos en este trabajo.
3.3 Simulación
Para comprender mejor el comportamiento de la varilla, hemos realizado una simulación en la que se representa una varilla de longitud [math] 1 m [/math] y se puede apreciar cómo se transmite el calor a lo largo del tiempo. Para ello, hemos elaborado el siguiente código de Matlab:
clear all
close all
% Coeficientes de Fourier
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
% Término k-ésimo de la solución del problema habiendo restado la solución estacionaria
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t));
% Definimos la solución estacionaria
sol_estacionaria=@(x) x;
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
% Definir el tiempo y el paso de tiempo
tiempo_final = 1;
paso_tiempo = 0.005;
tiempo = 0:paso_tiempo:tiempo_final;
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap1_5_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=14; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
title(['Varilla en t = ' num2str(t)]);
figura=figure(1);
% Calcular la temperatura en cada punto de la varilla en el tiempo t
w = zeros(size(Longitudinal));
for k = 1:1000
w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica/calor_especifico,c_k(Longitudinal(1,:),k));
end
u=w+sol_estacionaria(Longitudinal);
% Graficar la varilla con los colores según la temperatura
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title(['Varilla en t = ' num2str(t)]);
xlabel('Longitud');
axis([0 longitud -radio radio -radio radio]);
view(3);
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
end
% Cerrar el vídeo
close(pelicula);Al igual que en el otro código elaborado en esta sección, lo primero que definimos es la solución estacionaria del problema, así como los coeficientes de Fourier y los términos de la solución homogeneizada. Además, lo siguiente que debemos determinar son el coeficiente de conductividad térmica y el calor específico. En este caso, tomaremos ambos valores iguales a [math] 1 [/math]. Sin embargo, en la próxima sección estudiaremos el efecto que tienen estos valores en la solución.
A continuación, pasamos a introducir los datos que nos ayudarán a la representación de la varilla metálica, como son la longitud y el radio de esta o el número de puntos que pondremos a lo largo del eje x y el número de puntos que tomaremos para formar las circunferencias alrededor del eje longitudinal. Con estos datos generamos la malla de puntos de la varilla. Por otra parte, es fundamental definir los tiempos inicial y final, así como la longitud del intervalo temporal.
Para poder almacenar el resultado en un vídeo, hacemos uso de los comandos correspondientes propios de Matlab, entre los que destacan VideoWriter y FrameRate. Para cada uno de los tiempos definidos, debemos sumar los k-términos de la solución, siendo k el número que deseemos. A este resultado debemos añadirle la solución estacionaria para deshacer el cambio de variable y almacenamos el resultado en la variable u, que representa la temperatura de cada punto de la malla. Esta variable nos permitirá definir la escala de colores con la que vamos a representar la varilla metálica, por lo que a continuación, representamos la varilla haciendo uso de los datos de temperatura en cada punto para asignarle un color. Finalmente, añadimos la imagen generada al vídeo que se está creando y repetimos el proceso hasta haber recorrido todos los instantes de tiempo.
Cabe destacar que, el código realizado puede dar error en su ejecución si no se cuenta con una versión de Matlab de [math] 2022 [/math] o superior. En caso de no contar con dicha versión, basta sustituir en la línea 43 del código clim por caxis.
El resultado obtenido con este código es la siguiente simulación:
En esta simulación, apreciamos como la temperatura, partiendo de [math]0 °C [/math], aumenta de derecha a izquierda de acuerdo con las condiciones impuestas en el sistema de EDP’s y acorde a lo deducido del comportamiento del flujo en las gráficas anteriores.
4 Cambio de coeficientes
Planteamos ahora el mismo escenario anterior pero con coeficiente de conductividad térmica distintos: [math] k_1 =\frac{1}{2} \frac{cal}{ºC s m}, k_2=\frac{1}{6} \frac{cal}{ºC s m}[/math], en lugar de [math] k=1 \frac{cal}{ºC s m} [/math]. El calor específico [math] c [/math] continuará siendo igual a [math] 1 \frac{cal}{ºC g} [/math]. Por tanto, los sistemas quedarían como:
De manera análoga a los casos anteriores, obtenemos que la solución es de la forma
En primer lugar, es imprescindible recordar para la comparación que cuanto mayor sea el coeficiente de conductividad térmica [math] k [/math], mayor será la capacidad del material para conducir el calor. Además, estudiaremos cómo se comporta la solución del sistema homogeneizado en el punto intermedio de la barra, en [math] x = \frac{1}{2} [/math], y la compararemos para los coeficientes de conductividad térmica [math] k=1 [/math] y [math] k_1=\frac{1}{2} [/math]. Notar que, la solución del sistema homogeneizado equivale a [math] u \left(\frac{1}{2}, t\right) –v\left(\frac{1}{2}\right) [/math].
Para ello, hemos realizado el código presentado a continuación.
clear all
close all
% COEFICIENTES DE FOURIER
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t));
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) x;
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=[1,1/2,1/6];
calor_especifico=[1,1,1];
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.01;
punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=1;
paso_tiempo=0.0001;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X1,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
X1(:,end+1)=punto_espacio_estudio*ones(length(X1(:,1)),1); % vector con el punto del espacio a estudiar
X=transpose(sortrows(transpose(X1))); % añadimos el punto del espacio a estudiar en la matriz de mallado espacial
T(:,end+1)=T(:,end); % redimensionamos la matriz de mallado temporal
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
k_end=200; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
for coefi=1:length(coeficiente_conductividad_termica)
w_coef=zeros(size(X)); % inicializamos la variable w_coef1 (solución del problema con el primer coeficiente habiendo restado la solución estacionaria)
for ks=1:k_end % resolvemos el problema con los primeros k_end términos.
w_coef=w_coef+w_k(X,T,ks,coeficiente_conductividad_termica(coefi)/calor_especifico(coefi),c_k(X(1,:),ks));
end
u=w_coef+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
% REPRESENTACIONES DE LAS SOLUCIONES PARA LOS DOS COEFICIENTES DE CONDUCTIVIDAD
figure(1)
subplot(1,length(coeficiente_conductividad_termica),coefi)
mesh(X,T,u)
colormap('jet')
colorbar
xlabel('Espacio')
ylabel('Tiempo')
title("D = "+num2str(coeficiente_conductividad_termica(coefi)/calor_especifico(coefi)))
sgtitle("Representación de los primeros "+num2str(k_end)+" términos de la solución según los valores de D")
% REPRESENTAMOS LA DIFERENCIA ENTRE LA SOLUCIÓN Y LA SOLUCIÓN ESTACIONARIA EN punto_espacio_estudio
figure(2)
plot(T(:,1),w_coef(:,ceil(length(w_coef(1,:))/2)),'linewidth',1.5,'DisplayName',"Solución con D="+num2str(coeficiente_conductividad_termica(coefi)/calor_especifico(coefi)))
hold on
xlabel('Tiempo')
ylabel('Temperatura')
title("Diferencia entre la solución y la solución estacionaria en x="+num2str(punto_espacio_estudio)+" con respecto al tiempo")
end
figure(2)
legend('Location','southeast')
De nuevo, en primer lugar se definen los coeficientes de Fourier de la extensión impar de la función [math] x [/math] (solución estacionaria) como [math] c_k [/math], además del k-ésimo término de la solución homogeneizada y la solución estacionaria del sistema.
Posteriormente se incluyen distintos datos o parámetros modificables según las condiciones que deseemos establecer. Entre ellos se encuentran los valores de los coeficientes de conductividad térmica [math] k [/math], en forma de vector, que busquemos comparar.
Por otro lado, tenemos el dominio espaciotemporal que busquemos establecer. Al igual que en todos los códigos anteriores, de cara a aproximar los resultados continuos de una manera discreta, debemos establecer una malla espaciotemporal. Cuantos más puntos tenga esta malla, mayor precisión tendrá la aproximación al caso continuo.
Es destacable la definición de punto_espacio_estudio. Este punto es un punto específico dentro del dominio espacial en el cual queremos realizar el estudio. En dicho punto calcularemos cómo evoluciona la solución homogeneizada a medida que avanza el tiempo. Para añadir este punto estudio a la malla espaciotemporal, añadimos una columna con el punto repetido en todas las filas a dicha matriz de mallado, y posteriormente ordenamos la matriz por columnas.
Por último, realizamos una construcción análoga a la ya mencionada en el código de la sección anterior, a diferencia de que el proceso se debe repetir para cada uno de los coeficientes [math] k [/math] para los que deseemos conocer la solución.
Con respecto a esta primera imagen, como se puede observar, los dos primeros casos son muy similares, aunque hayamos cambiado el coeficiente de conductividad térmica ligeramente. Sin embargo, si comparamos los casos [math] k=1 \frac{cal}{ºC s m} [/math] y [math] k_2=\frac{1}{6} \frac{cal}{ºC s m}[/math], sí se observa una diferencia en el comportamiento. La superficie representada con [math] k_2= \frac{1}{6} \frac{cal}{ºC s m} [/math] presenta una mayor curvatura en el instante de tiempo final representado ([math] t=1 s[/math]).
Esto se debe a que aunque las tres soluciones se aproximan a la solución estacionaria [math] v(x,t)=x [/math], el coeficiente de conductividad térmica influye en la velocidad de convergencia de la solución hacia la solución estacionaria. Por ello, para un instante de tiempo [math] t=1 s[/math] se observan diferentes curvaturas en las superficies. En caso de aumentar el intervalo de tiempo lo suficiente, las diferencias en el instante final serían despreciables. Se anima al lector a cambiar el valor del parámetro tiempo_final por [math] 3’5 s[/math] en el código recién explicado.
Esta conclusión se puede observar perfectamente en esta segunda gráfica, en la que se muestra la solución del sistema homogeneizado en [math] x=\frac{1}{2} [/math] con los diferentes coeficientes de conductividad térmica. Como en la gráfica se representa la diferencia entre la solución y la estacionaria, cuanto más cercanos a cero sean estos valores, menos tiempo quedará para la convergencia de la solución.
Si prestamos atención a la gráfica, se puede observar a simple vista como, a medida que disminuimos los valores de los coeficientes de conductividad térmica [math] k [/math], el tiempo necesario para llegar a una temperatura nula es cada vez mayor. Luego, la velocidad de convergencia de la solución a la estacionaria aumenta a medida que aumenta el coeficiente de conductividad térmica.
Todo ello tiene sentido con el concepto de coeficiente de conductividad térmica dado que, los coeficientes de conductividad térmica altos conducen mejor el calor y, por tanto, tardan una menor cantidad de tiempo en estabilizar la temperatura del cuerpo con respecto a la temperatura ambiental.
4.1 Simulación
En esta sección, también se ha elaborado un código para la simulación y visualización del problema en el caso de la varilla metálica.
clear all
close all
% Coeficientes de Fourier
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
% Término k-ésimo de la solución del problema habiendo restado la solución estacionaria
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t));
% Definimos la solución estacionaria
sol_estacionaria=@(x) x;
% Valores de los coeficientes
coeficiente_conductividad_termica=[1,1/2,1/6];
calor_especifico=[1,1,1];
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
% Definir el tiempo y el paso de tiempo
tiempo_final = 1;
paso_tiempo = 0.005;
tiempo = 0:paso_tiempo:tiempo_final;
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap6_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=14; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
figura=figure(1);
sgtitle("Varilla en t = "+num2str(t)+" para distintos valores de D")
figura.Position=[100,100,1550,550]
for coefi=1:length(coeficiente_conductividad_termica)
% Calcular la temperatura en cada punto de la varilla en el tiempo t
w = zeros(size(Longitudinal));
for k = 1:1000
w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica(coefi)/calor_especifico(coefi),c_k(Longitudinal(1,:),k));
end
u=w+sol_estacionaria(Longitudinal);
% Graficar la varilla con los colores según la temperatura
subplot(1,length(coeficiente_conductividad_termica),coefi)
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title("D="+coeficiente_conductividad_termica(coefi)/calor_especifico(coefi));
xlabel('Longitud');
axis([0 longitud -radio radio -radio radio]);
view(3);
end
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
end
% Cerrar el vídeo
close(pelicula);
La construcción del código es análoga a la ya mencionada en la sección 3.3, pero con una pequeña diferencia. Se define un vector de coeficientes de conductividad térmica [math] k [/math] y se repite el proceso de representación para cada uno de los valores.
En este caso, se aprecian perfectamente las observaciones realizadas. El coeficiente de conductividad térmica afecta principalmente a la velocidad con la que se transmite el calor. Si observamos el primer y tercer video, la velocidad a la que se transmite el calor es muy distinta. En el primer caso, con coeficiente de difusión [math] D=1 [/math], el sistema alcanza el estado estacionario mucho más rápido. De hecho, el caso en el que el coeficiente de difusión [math] D=1/6 [/math], el sistema aún no ha alcanzado el estado estacionario cuando finaliza el vídeo.
5 Cambio de condiciones
Una vez hemos estudiado cómo influye el realizar un cambio en el coeficiente de conductividad térmica en el problema original, nos planteamos cómo varía la solución si suponemos que la temperatura en el extremo derecho es también de [math] 0 °C[/math] pero inicialmente la temperatura de la varilla viene dada por la función [math] u_0(x)=max\{0,1-4|x-1/2|\}[/math].
Por tanto, tendríamos el siguiente sistema:
La solución general de este sistema permanece igual que en sistema original a diferencia de que, en este caso, al tener las condiciones frontera homogeneizadas, la solución particular es nula. Es decir, la solución general es de la forma:
La única diferencia es la determinación de los coeficientes de Fourier [math] c_k [/math], en este caso, resulta mucho más sencillo aproximar el valor de las integrales de dichos coeficientes por la fórmula del trapecio.
A continuación, se presenta un código en Matlab que sigue la misma estructura que el presentado en las secciones [math] 3.2 [/math]. A diferencia de este, al haber cambiado la condición inicial por la función [math] u_0(x) [/math] y la condición frontera en el extremo derecho a una temperatura inicial nula, únicamente cambia que sol_estacionaria(x) es un vector de ceros y la determinación de los coeficientes de Fourier, en este caso calculados por la fórmula del trapecio.
clear all
close all
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(0,1-4.*abs(x-1/2*ones(size(x)))).*sin(k*pi.*x))./trapz(x,(sin(k*pi.*x)).^2);
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t));
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;
punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=1;
paso_tiempo=0.001;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
w=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
k_end=400; %último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
for ks=1:k_end % queremos resolver el problema para distintos valores de k
w=w+w_k(X,T,ks,coeficiente_conductividad_termica/calor_especifico,c_k(X(1,:),ks));
u=w+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
if ks <= 10 % representamos las gráficas de a lo sumo los 10 primeros términos de la solución
figure(1)
subplot(2,min(ceil(k_end/2),5),ks)
hold on
mesh(X,T,u)
colormap('jet')
colorbar
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
end
% COMPARACIÓN DE LA SOLUCIÓN PARA k=10 Y k=k_end
if (k_end>10) & (ks==10)
figure(2)
sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
subplot(1,2,1)
mesh(X,T,u)
colormap('jet')
colorbar
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
sgtitle('Representación del flujo con respecto al tiempo en los extremos para distintos valores de k')
subplot(1,2,1)
deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
plot(T(:,1),flujo0,'r','linewidth',1.5)
hold on
plot(T(:,1),flujo1,'b','linewidth',1.5)
xlabel('Tiempo')
ylabel('Flujo')
legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
title("k="+num2str(ks))
elseif (ks==k_end)
figure(2)
if (ks>10)
subplot(1,2,2)
end
mesh(X,T,u)
colormap('jet')
colorbar
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
if (ks>10)
subplot(1,2,2)
end
deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
plot(T(:,1),flujo0,'r','linewidth',1.5)
hold on
plot(T(:,1),flujo1,'b','linewidth',1.5)
xlabel('Tiempo')
ylabel('Flujo')
legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
title("k="+num2str(ks))
end
end
figure(4)
scatter(X(1,:),u(1,:),5, u(1,:),'filled')
colormap('jet')
title("Solución para k="+ num2str(k_end)+ " y t=0")
xlabel('Espacio')
ylabel('Temperatura')
Cómo se puede observar, la solución obtenida en la primera representación cumple con las condiciones establecidas en el sistema. Además, a medida que aumenta la precisión de la solución (introduciendo un mayor número de términos en la suma), esta alcanza un máximo de temperatura mayor y en un intervalo espacial más pequeño con respecto al punto medio de la varilla. Por otro lado, es destacable como la solución es estrictamente positiva si [math] t\gt0 [/math]. Esto se puede interpretar como que la difusión transporta energía con velocidad infinita de propagación.
La mejor manera de entender este fenómeno es observando esta segunda gráfica, en la que se ha representado la solución para los [math] k=10 [/math] y los [math] k=400 [/math] primeros términos. La solución para [math] k=400 [/math] tiene forma de pico. Esto es llamativo dado que la solución no es diferenciable en el instante inicial [math] t=0 [/math].
Sin embargo, en cuanto [math] t [/math] alcanza valores positivos, el calor que se encontraba acumulado en la zona central de la varilla se difunde a los extremos de la misma, provocando que la temperatura en los puntos en los que era nula en el instante inicial, pase a no serla. Cómo esto sucede para cualquier incremento de tiempo tan pequeño como deseemos, esto indica que la velocidad de propagación de calor es infinita.
En relación con la tercera gráfica obtenida con el código citado anteriormente, se puede observar como el flujo en ambos extremos comienza en 0. Esto se debe a que la condición inicial verifica que en los puntos próximos a dichos extremos la temperatura es también cero. Sin embargo, a medida que el calor concentrado alrededor del punto medio se dispersa hacia los extremos de la varilla, el flujo en ambos extremos se ve aumentado en valor absoluto. Cabe destacar que el flujo en el extremo izquierdo de la varilla (representado en color rojo) es negativo y, por ello, saliente; mientras que el flujo en el extremo derecho (en color azul) es positivo y, por ello, también saliente. El hecho de que ambos flujos en los extremos sean salientes indica que la varilla se enfría en su interior. Además, ambos flujos tienden a cero cuando se alcanza la solución estacionaria, lo que indica que la temperatura en la varilla es constante a lo largo de toda ella.
De nuevo, hemos implementado un código en Matlab de la varilla metálica para mejorar su comprensión. El código realizado es idéntico al realizado en la sección 3.3, únicamente se diferencian en las funciones y condiciones iniciales que las establecemos acordes a este nuevo problema. El código sería el siguiente:
clear all
close all
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(0,1-4.*abs(x-1/2*ones(size(x)))).*sin(k*pi.*x))./trapz(x,(sin(k*pi.*x)).^2);
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t));
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
% Definir el tiempo y el paso de tiempo
tiempo_final = 0.55;
paso_tiempo = 0.001;
tiempo = 0:paso_tiempo:tiempo_final;
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap7_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=37; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
figura=figure(1);
% Calcular la temperatura en cada punto de la varilla en el tiempo t
w = zeros(size(Longitudinal));
for k = 1:200
w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica/calor_especifico,c_k(Longitudinal(1,:),k));
end
u=w+sol_estacionaria(Longitudinal);
% Graficar la varilla con los colores según la temperatura
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title(['Varilla en t = ' num2str(t)]);
xlabel('Longitud');
axis([0 longitud -radio radio -radio radio]);
view(3);
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
if t==0
imwrite(imagen.cdata,"Ej1_Ap7_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
end
end
% Cerrar el vídeo
close(pelicula);
Como se puede observar, la varilla comienza con una temperatura inicial más alta en el interior, pero con el trascurso del tiempo, se va enfriando hasta alcanzar el estado estacionario con temperatura nula en toda la varilla.
Por último, dado que la condición inicial establecida (a la que hemos denotado como [math] u_0(x) [/math]), puede resultar algo difícil de visualizar, hemos realizado una representación gráfica de dicha varilla en el instante inicial. Con esta representación gráfica podemos hacernos a la idea de cómo se distribuye la temperatura en [math] t=0 [/math] de una manera más sencilla para así acabar de entender este caso.
6 Cambios de condiciones de frontera
Ahora, en vez de suponer que la temperatura en el extremo derecho es de [math]0ºC [/math], suponemos que la barra está aislada térmicamente en dicho extremo. El sistema de EDP’s quedaría de la siguiente manera:
Resolviendo el sistema por separación de variables obtenemos como solución la serie
[math] u(x,t)=\sum_{k=1}^\infty c_{k} sin\left(\left(k \pi-\frac{\pi}{2}\right)x\right)e^{(k \pi -\frac{\pi}{2})^2 Dt }+x [/math]. Esta solución no queda en función de la base trigonométrica usual, no obstante el conjunto de funciones [math]\left\{ sin\left(\left(k \pi-\frac{\pi}{2}\right)x\right)\right\} _{k \in \mathbb{N}} [/math] es un conjunto completo en [math] x\in [0,1] [/math] y es un conjunto ortogonal:
[math] \int^{1}_{0} \sin\left(\left(n \pi-\frac{\pi}{2}\right)x\right) \sin\left(\left(m \pi-\frac{\pi}{2}\right)x\right)dx= \int^{1}_{0} \frac{\cos\left(\left(n\pi-m\pi\right)x\right) -\cos\left(\left(n\pi+m\pi - \pi\right)x\right)}{2}=0[/math], para [math] n,m \in \mathbb{N}[/math]
Por consiguiente, podemos aproximar el dato inicial usando los coeficientes de Fourier asociados a esta nueva base. Calculando estos coeficientes en Matlab mediante la fórmula del trapecio, obtenemos el correspondiente comportamiento de la solución para distinto número de términos de la solución, es decir, valores de k:
clear all
close all
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(zeros(size(x)),1-4.*abs(x-1/2*ones(size(x)))).*sin((k*pi-pi/2).*x))./trapz(x,(sin((k*pi-pi/2).*x)).^2);
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin((k*pi-pi/2).*x).*exp(-D.*((-pi/2+k*pi)^2.*t));
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;
punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=2;
paso_tiempo=0.001;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
w=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
k_end=1000; %último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
for ks=1:k_end % queremos resolver el problema para distintos valores de k
w=w+w_k(X,T,ks,coeficiente_conductividad_termica/calor_especifico,c_k(X(1,:),ks));
u=w+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
if ks <= 10 % representamos las gráficas de a lo sumo los 10 primeros términos de la solución
figure(1)
subplot(2,min(ceil(k_end/2),5),ks)
hold on
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
end
% COMPARACIÓN DE LA SOLUCIÓN PARA k=10 Y k=k_end
if (k_end>10) & (ks==10)
figure(2)
sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
subplot(1,2,1)
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
sgtitle('Representación del flujo con respecto al tiempo en los extremos para distintos valores de k')
subplot(1,2,1)
deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
plot(T(:,1),flujo0,'r','linewidth',1.5)
hold on
plot(T(:,1),flujo1,'b','linewidth',1.5)
xlabel('Tiempo')
ylabel('Flujo')
legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
title("k="+num2str(ks))
elseif (ks==k_end)
figure(2)
if (ks>10)
subplot(1,2,2)
end
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
titulo="k = "+num2str(ks);
title(titulo)
figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
if (ks>10)
subplot(1,2,2)
end
deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
plot(T(:,1),flujo0,'r','linewidth',1.5)
hold on
plot(T(:,1),flujo1,'b','linewidth',1.5)
xlabel('Tiempo')
ylabel('Flujo')
legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
title("k="+num2str(ks))
end
end
En primer lugar, para valores suficientemente altos de k, en el espacio [math]x=0 [/math] la solución se mantiene constantemente nula como indica la condición frontera. Por otra parte, en los espacios intermedios de la varilla para [math] t=0 [/math] la solución toma valores más elevados acorde a la condición inicial [math] u(x,0)= max\{0,1-4|x-1/2|\}, \quad 0 \lt x \lt 1[/math]. Sin embargo, vamos a aumentar el número de términos hasta[math]k=1000[/math] para analizar con mayor precisión el fenómeno que sucede en el extremo derecho
En comparación con el caso anterior, el calor no se transmite en ambos extremos de la varilla de manera uniforme, sino que hay mayor transporte de calor hacia el lado izquierdo. Esto se debe a que hemos impuesto que el extremo derecho esté aislado térmicamente. Observando la gráfica con detalle, se aprecia que en el extremo izquierdo hay una clara diferencia de temperatura con respecto a los puntos próximos, por lo que en este extremo sí habrá flujo. Sin embargo, cerca del extremo derecho la temperatura se mantiene constante como consecuencia directa de que la derivada espacial en ese extremo sea nula, es decir, el flujo sea nulo. Para comprobar estas ideas, representamos el flujo en ambos extremos de la varilla:
Podemos ver al igual que en la otra gráfica como en el extremo derecho no existe flujo y en el extremo izquierdo, el flujo será saliente. Como vimos en la primera gráfica, el flujo empieza en [math] 0 [/math] debido a su condición inicial. Rápidamente aumenta a valores más negativos acorde al aumento de la pendiente. En tiempo [math] t=0.2 ~s[/math] el flujo se reduce rápidamente, hasta que en [math] t=2~s [/math] el flujo desaparece y se llega a la solución estacionaria.
6.1 Simulación
Como hemos hecho en los casos anteriores, modelizamos el cambio de temperatura que se produce en la varilla:
clear all
close all
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(zeros(size(x)),1-4.*abs(x-1/2*ones(size(x)))).*sin((k*pi-pi/2).*x))./trapz(x,(sin((k*pi-pi/2).*x)).^2);
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin((k*pi-pi/2).*x).*exp(-D.*((-pi/2+k*pi)^2.*t));
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
% Definir el tiempo y el paso de tiempo
tiempo_final = 2;
funcion_velocidad_t=@(t) t*1.05;
tiempo=[0,0.001];
while tiempo(end)<tiempo_final
tiempo(end+1)=funcion_velocidad_t(tiempo(end));
end
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap8_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=11; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
figura=figure(1);
% Calcular la temperatura en cada punto de la varilla en el tiempo t
w = zeros(size(Longitudinal));
for k = 1:1000
w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica/calor_especifico,c_k(Longitudinal(1,:),k));
end
u=w+sol_estacionaria(Longitudinal);
% Graficar la varilla con los colores según la temperatura
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title(['Varilla en t = ' num2str(t)]);
xlabel('Longitud');
axis([0 longitud -radio radio -radio radio]);
view(3);
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
if t==0
imwrite(imagen.cdata,"Ej1_Ap8_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
end
end
% Cerrar el vídeo
close(pelicula);
Esta animación nos afirma todas las conclusiones obtenidas a raíz de las otras gráficas. En concreto, se puede apreciar con claridad como cerca del extremo derecho la temperatura se mantiene constante. Con respecto a los cambios realizados en este código con respecto a las anteriores simulaciones, se han modificado los tiempos para que el video capture más imágenes al principio, donde ocurren cambios significativos con mayor frecuencia, y luego tome imágenes en intervalos más espaciados de tiempo, ya que los cambios son menos notorios.
6.2 Principio del máximo
Por último, vamos a comprobar que el principio del máximo se verifica para este problema, con las condiciones tomadas en el apartado anterior.
Para este último caso, el espacio es [math]Q_{T}=(0,1)\times(0,2) [/math] y se verifica que la solución [math] u(x,t)=\sum_{k=1}^\infty c_{k} sin\left(\left(k \pi-\frac{\pi}{2}\right)x\right)e^{(k \pi -\frac{\pi}{2})^2 Dt }+x :=\sum_{k=1}^\infty u_k (x,t)[/math] converge uniformemente y las derivadas de [math]u_k (x,t)[/math] también. Por tanto, podemos afirmar que:
[math] u_x(x,t) =\sum_{k=1}^\infty \frac{d}{dx} u_k (x,t) [/math]
[math] u_t(x,t) =\sum_{k=1}^\infty \frac{d}{dt}u_k (x,t) [/math]
Como podemos derivar para cualquier orden, se concluye que [math] u(x,t) \in C^{\infty}(Q_{T}) [/math].
Además, [math] u_t(x,t)-u_{xx}(x,t)= \sum_{k=1}^\infty \frac{d}{dt} u_k (x,t)- \frac{d^2}{dx^2}u_k (x,t) =0 \leq 0[/math]. Como las condiciones del principio del máximo se verifican, el máximo de la función ha de estar en la frontera parábolica de [math] Q_T [/math]. En nuestro caso, se puede apreciar como el máximo de la solución es el valor [math]1[/math], que se alcanza en el tiempo inicial[math]t=0 s[/math] y espacio [math]x=0.5 m[/math], perteneciendo claramente a la frontera parábolica de nuestro problema.
7 Solución fundamental de la ecuación del calor
En esta sección vamos a estudiar una solución al sistema
Para ello, vamos a usar la solución fundamental de la ecuación del calor. Cabe destacar que como este sistema no tiene condiciones frontera, pues [math] x \in \mathbb{R} [/math], es posible que no haya unicidad de las soluciones.
7.1 Definición solución fundamental
La solución fundamental de la ecuación del calor viene dada por
Si el lector lo desea, puede consultar la deducción de esta solución en la página 42 del libro Partial differential equations in action from modelling to theory de Sandro Salsa.
Destacamos que, buscando la conservación de la masa total para cualquier tiempo, es lógico que se verifique que [math]\int_{\mathbb{R^{n}}} u(x,t) d\mathbf{x} =q[/math] , con [math] q [/math] constante para todo tiempo. Sin embargo, en la deducción de esta fórmula se ha impuesto la condición [math] q=1 [/math], dando como resultado una familia de funciones Gaussianas parametrizadas en el tiempo.
7.2 Caso unidimensional
En concreto, en el caso unidimensional, la solución fundamental de la ecuación del calor es
Tomando el valor de la constante de difusión [math] D=1 [/math], [math] x \in [-1,1] [/math] y [math] t \in [10^{-2},1] [/math]; hemos elaborado el siguiente código de Matlab para representar la solución
clear all
close all
%DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
%VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
conduccion_calor=1;
%DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-1;
extremo_derecho=1;
paso_espacio=0.001;
%DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=10^(-2);
tiempo_final=1;
paso_tiempo=0.001;
%DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
%REPRESENTAMOS LA SOLUCIÓN FUNDAMENTAL
mesh(X,T,sol_fundamental(X,T,coeficiente_conductividad_termica/conduccion_calor))
En este código, lo primero que definimos es la solución fundamental en una dimensión. A continuación, simplemente añadimos los valores del coeficiente de conductividad térmica y del calor específico, que en este caso ambos son iguales a [math] 1 [/math]. Además, al igual que sucedía en los códigos anteriormente explicados, en Matlab trabajamos de manera discreta. Por ello, lo siguiente a definir son los datos que vamos a necesitar para la discretización del espacio y del tiempo y el mallado espaciotemporal haciendo uso de la función de Matlab meshgrid, para así poder hacer la representación gráfica de la solución con respecto al espacio y al tiempo con la función de Matlab mesh.
En la gráfica obtenida como resultado, se aprecian claramente algunas propiedades. La primera de ellas es la paridad de la función en la variable [math] x [/math]. De hecho, de manera analítica se puede comprobar que:
Demostrándose así la simetría. Además, lo más destacable de la gráfica es el fenómeno que se produce cuando nos aproximamos a [math] t=0 [/math]. Para poder realizar un mejor estudio de este fenómeno, hemos representado la solución fundamental para [math] x \in [-10,10] [/math] y para [math] t \in [10^{-2},5] [/math].
En esta gráfica, se aprecia perfectamente que cuando [math] t [/math] tiende a [math] 0 [/math], la solución fundamental tiende a [math] 0 [/math] en cualquier [math] x \neq 0 [/math] y tiende a [math] \infty [/math] en [math] x=0 [/math]. Esta distribución se conoce como la Delta de Dirac.
Se define la Delta de Dirac como el límite de una sucesión de funciones [math] \{f^{n}\}_{n \in \mathbb{N}} [/math] que tienden a cero en todo punto del espacio excepto en un punto para el cual diverge a infinito. Concretamente, llamaremos Delta de Dirac al siguiente límite: [math] \delta = \lim_{n \to \infty} f^n(x) [/math] con [math] \lim_{n \to \infty} f^n(0)=\infty [/math] y [math] \lim _{n \to \infty} f^n(x)=0 [/math] para [math] x \neq 0[/math] que además verifica [math] \int_{\mathbb{R}} f^n(x) \cdot dx =1 [/math].
Además, vemos como el máximo de la solución fundamental se alcanza en [math] (x,t)=(0,0) [/math], cumpliéndose así el principio del máximo. Además, el valor de la función disminuye rápidamente a medida que nos alejamos de este punto en ambas direcciones. Finalmente, cabe destacar que la solución esatacionaria correspondiente a este problema es [math] v(x)=0 [/math]. De hecho,
Sin embargo, la solución fundamental nunca se anula si [math] t\gt0 [/math]. Esto quiere que el sistema tiende a la solución estacionaria, pero no la alcanza en tiempo finito.
7.2.1 Simulación
A continuación, se presenta un código correspondiente a la simulación de una varilla metálica sobre la que se difunde el calor según la solución fundamental. El funcionamiento y construcción de este código es análogo a los ya mencionados en otras secciones. El código sería el siguiente:
clear all
close all
% Definimos la solución fundamental
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
% Definir la longitud y radio de la varilla
longitud_positiva = 1;
radio = 0.25;
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
% Definir el tiempo y el paso de tiempo
tiempo_final = 2;
paso_tiempo = 0.005;
tiempo = 10^(-3):paso_tiempo:tiempo_final;
% Generar la malla de puntos en la varilla
longitudinal = linspace(-longitud_positiva, longitud_positiva, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej2_Ap1_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=27; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
figura=figure(1);
% Calcular la temperatura en cada punto de la varilla en el tiempo t
u=sol_fundamental(Longitudinal,t,coeficiente_conductividad_termica/calor_especifico);
% Graficar la varilla con los colores según la temperatura
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title(['Varilla con temperatura en t = ' num2str(t)]);
xlabel('Espacio');
axis([-longitud_positiva longitud_positiva -radio radio -radio radio]);
view(3);
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
if t==0
imwrite(imagen.cdata,"Ej2_Ap1_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
end
end
% Cerrar el vídeo
close(pelicula);
Este caso se da en la vida real cuando el radio de la varilla es despreciable en comparación con la longitud de esta. Sin embargo, en esta simulación, hemos representado una varilla de longitud [math] 2m [/math] para que coincidan los resultados con la primera gráfica de la solución fundamental expuesta.
Como se puede observar, el tiempo inicial del video es [math] t=10^{-2} ~s[/math], puesto que la solución fundamental tiende a la delta de Dirac cuando [math] t \rightarrow 0 [/math], y por lo tanto, la varilla tendría temperatura infinita en sus puntos centrales y cero en el resto de puntos. Por otro lado, se puede ver como las temperaturas más elevadas inicialmente se distribuyen en el centro de la varilla y con el tiempo se van difundiendo hacia los extremos tendiendo la temperatura de la varilla a [math] 0 ºC [/math], aunque nunca alcanzando el estado estacionario.
7.3 Resolución de la ecuación del calor con la solución fundamental
En esta sección, buscamos resolver el siguiente problema:
Llamamos soluciones autosemejantes a aquellas para las que existe [math] \alpha \in \mathbb{R} [/math] tal que [math] u(\lambda x, \lambda^2 t)= \lambda ^{\alpha} u(x,t) \quad \forall \lambda \gt 0, \forall (x,t) \in Q_{T} [/math].
En primer lugar, vamos a buscar las soluciones autosemejantes de la ecuación del calor [math] u_t-u_{xx}=0 [/math] de la forma [math] u(x,t)=U\left( \frac{x}{\sqrt{t}} \right) [/math]. De esta manera, se tiene que
Y, por lo tanto,
Llamando [math] \xi=\frac{x}{\sqrt{t}} [/math], podemos simplificar al expresión
Definiendo [math] f(\xi)= U’\left( \xi \right) [/math] y aplicando el método de separación de variables, se llega a que
A continuación, aplicando el teorema fundamental del cálculo, se deduce que
Para calcular los valores de [math] D,E \in \mathbb{R} [/math], debemos hacer uso de las condiciones iniciales y de frontera de las que disponemos. En primer lugar, se debe cumplir que [math] u(0,t)=1 [/math] o, lo que es lo mismo,
Por lo tanto, [math] U(\xi)=D \int_{0}^{\xi} e^{-\frac{\eta^2}{4}} d\eta +1, \quad D \in \mathbb{R} [/math]. Por otra parte, el ejercicio exige que [math] u(x,t) \rightarrow 0 [/math] cuando [math] x \rightarrow \infty [/math]. Esto es equivalente a que [math] U(\xi) \rightarrow 0 [/math] cuando [math]\xi \rightarrow \infty [/math] y haciendo uso de la paridad de la función [math] e^{-\frac{\eta^2}{4}}, ~\eta \in \mathbb{R} [/math], se obtiene la siguiente serie de igualdades:
Además, si aplicamos el cambio de variable [math] z=\frac{\eta}{2} [/math] y sabiendo que [math] \int_{\mathbb{R}} e^{-z^2} dz= \sqrt{\pi} [/math], llegamos a que
De esta manera, concluimos que la solución buscada es
siendo [math] erf(x)= \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-z^2} dz [/math].
Por lo tanto, deshaciendo el cambio de variable,
Cabe destacar que esta solución también verifica la condición inicial que faltaba por verificar, pues [math] \lim_{t \to 0} u(x,t)= \lim_{t \to 0} 1- \frac{2}{\sqrt{\pi}}\int_{0}^{\frac{x}{2 \sqrt{t}}} e^{-z^2} dz=1- \frac{2}{\sqrt{\pi}}\frac{1}{2}\int_{\mathbb{R}} e^{-z^2} dz=1- \frac{2}{\sqrt{\pi}}\frac{1}{2} \sqrt{\pi}=0 [/math].
Además,
Sin embargo, de nuevo ocurre que la función solución nunca toma el valor 1, pues [math] erf\left(\frac{x}{2\sqrt{t}}\right) [/math] nunca se anula. Esto indica que a pesar de que el sistema tienda a la solución estacionaria <v(x)=1>, esta no se alcanza en tiempo finito.
Una vez calculada la solución al problema, hemos elaborado un código en Matlab para representarla:
clear all
close all
% SOLUCIÓN DEL EJERCICIO 10b
sol_erfc=@(x,t) 1-erf(x./(2.*sqrt(t)));
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=20;
paso_espacio=0.1;
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=10;
paso_tiempo=0.001;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
% REPRESENTAMOS LA SOLUCIÓN
figure(1)
mesh(X,T,sol_erfc(X,T))
title('Representación de la solución')
xlabel('Espacio')
ylabel('Tiempo')
zlabel('Temperatura')
En este código, primero hemos definido la solución del ejercicio como una función. A continuación, como en el resto de los programas, hemos definido los datos necesarios para poder hacer una discretización del tiempo y del espacio. En este caso, se ha representado la función en [math] (x,t) \in [0,20] \times [0,10] [/math]. Con respecto a la gráfica obtenida como resultado, se puede apreciar a la perfección que se cumple que [math] u(0,t)=1, t\gt0 [/math] y [math] u(x,t) \rightarrow 0 [/math] cuando [math] x \rightarrow \infty [/math]. Además, [math] u(x,0)=0, x\gt0 [/math]. En esta gráfica, es importante destacar que se produce una discontinuidad en [math] (x,t)=(0,0)[/math]. Esta discontinuidad se debe a que las propias condiciones iniciales y de frontera no coinciden en dicho punto. A pesar de esto, Matlab hace una representación continua de la función.
7.3.1 Simulación
De nuevo, hemos hecho la simulación en una varilla metálica para poder ver visualmente cómo se comporta el calor. Como en este caso la solución está definida para todo [math] \forall x\gt0 [/math], hemos representado un trozo de la varilla metálica con una longitud de [math] 100m [/math] con un radio [math] 0.25m [/math].
clear all
close all
% SOLUCIÓN DEL EJERCICIO 10b
sol=@(x,t) 1-erf(x./(2.*sqrt(t)));
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
% Definir la longitud y radio de la varilla
longitud = 100;
radio = 0.25;
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 5001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
% Definir el tiempo y el paso de tiempo
tiempo_final = 500;
paso_tiempo = 1;
tiempo = 0:paso_tiempo:tiempo_final;
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej2_Ap2_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=34; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
figura=figure(1);
% Calcular la temperatura en cada punto de la varilla en el tiempo t
u=sol(Longitudinal,t);
% Graficar la varilla con los colores según la temperatura
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title(['Varilla con temperatura en t = ' num2str(t)]);
xlabel('Espacio');
axis([0 longitud -radio radio -radio radio]);
view(3);
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
if t==0
imwrite(imagen.cdata,"Ej2_Ap2_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
end
end
% Cerrar el vídeo
close(pelicula);
Este caso es muy distinto a los estudiados con anterioridad. El calor se concentra en el extremo izquierdo de la varilla en el instante inicial y con el paso del tiempo, se desplaza hacia la derecha constantemente. Recordamos que el radio de la varilla es despreciable en comparación con su longitud. Si tenemos en cuenta esto, podemos considerar dicha longitud infinita y, por ello, la solución estacionaria no se alcanza en tiempo finito. De hecho, destacamos que el vídeo recién introducido alcanza el tiempo [math] t=500~s[/math] y en dicho instante el calor se ha difundido muy poco en comparación con la longitud de la varilla.
7.4 Solución de la ecuación del calor empleando la convolución
Planteamos ahora un nuevo escenario. En este caso, partimos de la ecuación del calor [math] u_t-u_{xx}=0 [/math] en [math] x \in \mathbb{R} [/math] asociada al dato inicial [math] u_0(x)=1_{[-1,1]} [/math].
En este caso, la solución de dicha ecuación viene dada por la convolución:
Como la condición inicial es [math] u_0(x)=1_{[-1,1]} [/math], la integral de la convolución de la solución fundamental y la condición inicial será nula fuera del dominio de integración [math] [-1,1] [/math], por lo que es más sencillo reducir el cálculo de dicha convolución en todo [math] \mathbb{R} [/math] a únicamente el intervalo [math] [-1,1] [/math].
De esta manera, podemos dibujar la solución [math] u(x,t) [/math] en diferentes instantes de tiempo. Para ello, hemos implementado el siguiente código en Matlab:
clear all
close all
% DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
% DEFINIMOS LA FUNCIÓN SOLUCIÓN EN EL INSTANTE INICIAL
u_0=@(x) 1*((-1<=x)&(x<=1));
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% DEFINIMOS Y COMO LA VARIABLE SOBRE LA QUE VAMOS A INTEGRAR
Y=linspace(-1,1,100);
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-5;
extremo_derecho=5;
paso_espacio=0.001;
X=[extremo_izquierdo:paso_espacio:extremo_derecho];
% TIEMPOS A ESTUDIAR
T=[0.001,0.01,0.1];
t_final=100000;
% DEFINIMOS LA SOLUCIÓN u(x,t)
figure(1)
u=zeros(length(T),length(X));
for j=1:length(T)
for i=1:length(X)
u(j,i)=trapz(Y,sol_fundamental(X(i)-Y,T(j),coeficiente_conductividad_termica/calor_especifico).*u_0(Y));
end
subplot(1,length(T),j)
plot(X,u(j,:),'b','LineWidth',1.5)
axis([-2,2,0,1.001])
xlabel('Espacio')
ylabel('Temperatura')
title("t="+num2str(T(j)))
end
sgtitle('Representación de la solución para distintos instantes de tiempo t')
% Crear un vídeo
pelicula=VideoWriter('Solucion_reales','MPEG-4'); %creo el video
pelicula.FrameRate=15; %controla velocidad
open(pelicula); %lo abro
funcion_velocidad_t=@(t) t*1.1;
vector_t=[0,0.001];
while vector_t(end)<t_final
vector_t(end+1)=funcion_velocidad_t(vector_t(end));
end
for j=1:length(vector_t)
figura=figure(2);
for i=1:length(X)
u(j,i)=trapz(Y,sol_fundamental(X(i)-Y,vector_t(j),coeficiente_conductividad_termica/calor_especifico).*u_0(Y));
end
plot(X,u(j,:),'b','LineWidth',1.5)
axis([extremo_izquierdo,extremo_derecho,0,1.001])
view(2);
xlabel('Espacio')
ylabel('Temperatura')
title("t="+num2str(vector_t(j)))
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
end
close(pelicula);
Primero de todo, se define la solución fundamental y la condición inicial, es decir, la función [math] u_0(x) [/math] y los valores de los coeficientes de conductividad térmica y de calor específico.
Se define la variable [math] Y [/math] como la variable sobre la que vamos a integrar. A continuación, se discretiza el espacio definiendo por tanto los extremos donde está definida la solución y el paso de discretización deseado entre dichos puntos. Por último, antes de definir la solución se establece un vector de tiempos [math] T [/math] donde se quiere realizar un estudio de la solución obtenida en dichos instantes de tiempo.
Seguidamente, se define la solución dada por la convolución de las dos funciones ya mencionadas y se aproxima haciendo uso de la fórmula del trapecio en el intervalo espacial definido. En último lugar, se definen los comandos adecuados requeridos para la realización de un video que muestre las soluciones de la ecuación en distintos instantes de tiempo. Concretamente, desde [math]t= 0.0011 ~s [/math] a [math] t=97451 ~s[/math].
En relación con los resultados obtenidos al implementar el código, son destacables varios aspectos.
En primer lugar, como era de esperar, para valores de [math] t [/math] muy cercanos al valor [math] t=0 [/math] la solución [math] u(x,t) [/math] obtenida es muy similar a la condición inicial [math] u_0(x) [/math].
En segundo lugar, se puede apreciar como a medida que aumenta el tiempo, la solución [math] u(x,t) [/math] tiende a cero, la solución estacionaria de la ecuación, aunque no la alcanza en tiempo finito. Conceptualmente esto quiere decir que para tiempos lo suficientemente grandes se produce una variación despreciable de la temperatura en el objeto de estudio.
7.4.1 Simulación
Calculada la nueva solución, se presenta el código correspondiente a la simulación de la difusión del calor sobre la varilla. Su construcción es análoga a los anteriores.
clear all
close all
% DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
% DEFINIMOS LA FUNCIÓN SOLUCIÓN EN EL INSTANTE INICIAL
u_0=@(x) 1*((-1<=x)&(x<=1));
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
% DEFINIMOS Y COMO LA VARIABLE SOBRE LA QUE VAMOS A INTEGRAR
Y=linspace(-1,1,100);
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-longitud;
extremo_derecho=longitud;
paso_espacio=0.001;
X=[extremo_izquierdo:paso_espacio:extremo_derecho];
% Definir el tiempo y el paso de tiempo
tiempo_final = 10;
funcion_velocidad_t=@(t) t*1.05;
tiempo=[0.001];
while tiempo(end)<tiempo_final
tiempo(end+1)=funcion_velocidad_t(tiempo(end));
end
% DEFINIMOS Y COMO LA VARIABLE SOBRE LA QUE VAMOS A INTEGRAR
Y=linspace(-1,1,100);
% Generar la malla de puntos en la varilla
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(X, circunferencia);
% Crear un vídeo
pelicula=VideoWriter('Ej2_Ap3_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=13; %controla velocidad
open(pelicula); %abrimos el vídeo
% Bucle para cada paso de tiempo
for t = tiempo
figura=figure(1);
u=zeros(size(Longitudinal));
% Calcular la temperatura en cada punto de la varilla en el tiempo t
for i=1:length(Longitudinal(1,:))
u(1,i)=trapz(Y,sol_fundamental(X(i)-Y,t,coeficiente_conductividad_termica/calor_especifico).*u_0(Y));
end
for j=1:length(Longitudinal(:,1))
u(j,:)=u(1,:);
end
% Graficar la varilla con los colores según la temperatura
colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
colorbar;
clim([0, 1]); % Fijar el rango de colores
% Añadir etiquetas y ajustes de la figura
title(['Varilla con temperatura en t = ' num2str(t)]);
xlabel('Espacio');
axis([-longitud longitud -radio radio -radio radio]);
view(3);
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
if t==0
imwrite(imagen.cdata,"Ej2_Ap3_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
end
end
% Cerrar el vídeo
close(pelicula);
Como se puede observar en el vídeo, la varilla metálica está definida en el dominio espacial [math] [-1,1] [/math]. En el instante [math] t=0 [/math] se verifica la condición inicial y, por ello, toda la varilla comienza con una temperatura de [math] 1ºC [/math]. Con el transcurso del tiempo, el calor comienza a difundirse por el resto de la varilla. Es destacable que se ha establecido el tiempo de ejecución de manera que se reproduce mucho más lento y detallado los instantes de tiempo iniciales y, de una manera más rápida a medida que pasa el tiempo. Por último, la temperatura de la varilla tiende a estabilizarse a una temperatura nula.
7.5 Solución fundamental de la ecuación del calor en dimensión 2
Una vez hemos estudiado y analizado cómo se comporta la solución fundamental del calor en una dimensión vamos a analizar brevemente el caso de la solución fundamental en dimensión [math] 2 [/math]. La solución, tal y como hemos mencionado previamente, sería:
En dos dimensiones, la geometría del problema puede ser más representativa, especialmente en sistemas bidimensionales como placas o láminas. Por otro lado, el sistema y soluciones pueden resultar más complejos. Este aspecto gana importancia cuando tenemos casos con propiedades significativas del medio o condiciones de contorno o iniciales más elaboradas.
Asimismo, en dimensión [math] 2 [/math] es posible modelar fenómenos que no se pueden capturar en una dimensión, como la difusión del calor a lo largo de un plano sin necesidad de que el flujo sea uniforme en todos los puntos de una misma sección y que la difusión se produzca únicamente en una dirección fijada.
Por último, cabe destacar que la solución estacionaria de la solución fundamental es:
Sin embargo, de nuevo vuelve a ocurrir que la solución fundamental no se anula en ningún punto para t>0. Esto quiere decir que la solución fundamental tiende a la solución estacionaria, pero no la alcanza en tiempo finito.
La manera más sencilla de entenderlo es con una representación gráfica. Por tanto, hemos realizado un código en Matlab que implementa dicha solución en distintos instantes de tiempo. El código es el siguiente:
clear all
close all
% DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x1,x2,t,D) 1./((4.*pi.*D.*t)).*exp((-sqrt(x1.^2+x2.^2).^2)./(4.*D.*t));
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-1;
extremo_derecho=1;
paso_espacio=0.001;
[X1,X2]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,extremo_izquierdo:paso_espacio:extremo_derecho);
% TIEMPOS A ESTUDIAR
T=[0.00001,0.001,0.01,0.1];
t_final=0.001;
% DEFINIMOS LA SOLUCIÓN u(x,t)
figure(1)
for j=1:length(T)
subplot(1,length(T),j)
mesh(X1,X2,sol_fundamental(X1,X2,T(j),coeficiente_conductividad_termica/calor_especifico))
xlabel('Valores de x_1')
ylabel('Valores de x_2')
zlabel('Temperatura')
title("t="+num2str(T(j)))
end
sgtitle('Representación de la solución fundamental bidimensional para distintos instantes de tiempo t')
% Crear un vídeo
pelicula=VideoWriter('Solucion_fundamental_bidimensional','MPEG-4'); %creo el video
pelicula.FrameRate=10; %controla velocidad
open(pelicula); %lo abro
funcion_velocidad_t=@(t) t*0.95;
vector_t=[0.05];
while vector_t(end)>t_final
vector_t(end+1)=funcion_velocidad_t(vector_t(end));
end
for j=1:length(vector_t)
figura=figure(2);
mesh(X1,X2,sol_fundamental(X1,X2,vector_t(j),coeficiente_conductividad_termica/calor_especifico))
axis([-1,1,-1,1,0,90])
view(3);
xlabel('Valores de x_1')
ylabel('Valores de x_2')
zlabel('Temperatura')
colormap(jet);
colorbar
clim([0,90])
title("t="+num2str(vector_t(j)))
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
end
close(pelicula);
La construcción de este código es análoga a la ya realizada en la sección 7.2. La única diferencia sería la definición de la solución fundamental, ya que habría que introducir la correspondiente a dimensión [math] 2 [/math], además de ajustar la definición de los mallados, que ahora sería un mallada únicamente espacial. Además, se hace una gráfica en 3 dimensiones para cada valor del tiempo.
Como conclusión de estas representaciones se puede observar que cuanto más cerca es el valor de la [math] t [/math] a cero, más singular es la solución. Este fenómeno se aprecia a la perfección en el siguiente vídeo, en el que el tiempo se va aproximando a [math] 0 ~s [/math]:
Concretamente, al igual que en el caso unidimensional, la solución fundamental de la ecuación del calor en dimensión [math] n=2 [/math] coincide con la delta de Dirac (bidimensional) para [math] t=0 [/math]. Esta observación se aprecia perfectamente en el vídeo realizado. La delta de Dirac en dos dimensiones se define de manera similar a la delta de Dirac unidimensional:
Por otra parte, recordando que en el vídeo el tiempo decrece, se aprecia claramente cómo tiende a la solución estacionaria, [math] v(\mathbf{x})=0 [/math] cuando [math] t [/math] crece.
8 Referencias
- SectionB. Similarity solutions.
- Ecuaciones en derivadas parciales Facultad de Matemáticas. Universidad de Sevilla
- Mauro Milonga Miguel. Trabajo de fin de grado: Ecuación del calor. Universidad de Murcia
- [Partial differential equations in action from modelling to theory. Sandro Salsa]
