Ecuación del calor (Grupo GIXP)

Trabajo realizado por estudiantes
Título	Ecuación del calor. Grupo GIXP
Asignatura	EDP
Curso	2024-25
Autores	Gonzalo Garelly Israel López Francisco Lavao Paula León
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción Y enfoque

La modelización matemática es una herramienta fundamental en el ámbito de la física para predecir cómo actuará un sistema a lo largo del tiempo y el espacio y así poder actuar conforme a él o ajustarlo si lo permite. Un modelo muy conocido es el de la ecuación de difusión, que forma parte de las tres ecuaciones en derivadas parciales más empleadas en la física.

Dicha ecuación aparece en fenómenos como los siguientes:

• Gestión térmica en electrónica: permite modelar la disipación de calor en procesadores y optimizar el diseño de disipadores metálicos.
• Climatología: estudia la transferencia de calor en la atmósfera y los océanos, ayudando a predecir cambios climáticos.
• Biología: describe la difusión de fármacos en el torrente sanguíneo, facilitando el desarrollo de tratamientos más efectivos.

En este caso, se planteará el problema de la difusión del calor unidimensional en una varilla, para unas condiciones genéricas de coeficiente de difusión y de la función de generación de calor por fuerzas externas. Se modelizarán distintas situaciones, dadas unas condiciones de tipo Dirichlet e inicial fijadas, como usar varillas de distintos tipos de metal (variando el coeficientes de difusión) o aplicar calor en distintos puntos de la varilla mediante el contacto con hierros a distintas temperaturas.

2 Planteamiento del sistema de EDP

Imaginemos que tenemos una varilla metálica con densidad constante \( \rho \), de longitud \( L \), cuya temperatura inicial viene dada por la constante \( u_0 \). Además, los extremos de la varilla mantienen temperaturas fijas, con el extremo izquierdo siempre a una temperatura \( u_1 \) y el derecho a una temperatura \( u_2 \) (Condiciones Dirichlet).

También consideramos una función \( \omega \), que representa la producción de energía externa. Supondremos que es constante respecto del tiempo, es decir, \( \omega(x) \) con \( x \in [0,L] \).

Aplicando la ley de conservación de la energía, la ley de Fourier y el teorema de la divergencia, llegamos a la siguiente ecuación en derivadas parciales:

[math] u_t(x,t) - \frac{k}{Q \rho} u_{xx}(x,t) = \frac{1}{Q} \omega(x) [/math]

donde: \( k \) es la constante de conductividad térmica asociada a la ley de Fourier. \( Q \) es la capacidad calorífica, que relaciona la temperatura con la energía.

Reagrupando términos, definimos: [math] \alpha = \frac{k}{Q \rho}, \quad f(x) = \frac{1}{Q} \omega(x) [/math] donde \( \alpha \) se conoce como **coeficiente de difusión térmica**.

Con esto, si aplicamos nuestras condiciones en la frontera parabólica, el sistema queda planteado de la siguiente forma:

[math] \begin{cases} u_t - \alpha u_{xx} = f, & x \in [0,L], \quad t \gt 0 \\ u(x,0) = u_0, & x \in [0,L] \\ u(0,t) = u_1, & t \gt 0 \\ u(L,t) = u_2, & t \gt 0 \end{cases} [/math]

3 Resolución del sistema

Primero, buscamos la solución estacionaria \( v(x) \), que cumple:

[math] v_{xx} = -\frac{1}{\alpha} f(x) [/math]

Integrando dos veces:

[math] v(x) = -\frac{1}{\alpha} \int \int f(x) \, dx^2 [/math]

Imponiendo las condiciones de contorno:

[math] \begin{cases} v(0) = u_1 \\ v(L) = u_2 \end{cases} [/math]

Definimos \( w(x,t) = u(x,t) - v(x) \), lo que nos lleva a la ecuación homogénea:

[math] w_t - \alpha w_{xx} = 0 [/math] con condiciones frontera también homogéneas.

Aplicamos separación de variables y obtenemos la solución:

[math] w(x,t) = \sum_{n=1}^{\infty} c_n \sin\left(\frac{\sqrt{\alpha}n\pi}{L}x \right) e^{- \frac{\alpha n^2 \pi^2}{L^2} t} [/math]

Donde los coeficientes \( c_n \) los obtenemos al imponer que:

[math] w(x,0) = u_0 - v(x) = \sum_{n=1}^{\infty} c_n \sin\left(\frac{\sqrt{\alpha} n\pi}{L} x \right) [/math]

Realizando la extensión impar de \( u_0 - v(x) \) al intervalo \([-L, L]\) y hallando sus coeficientes en la serie de Fourier:

[math] c_n = \int_{-L}^{L} (u_0 - v(x)) \sin\left(\frac{\sqrt{\alpha} n\pi}{L} x \right) dx [/math]

Finalmente nos acordamos de que u(x,t)=w(x,t)+v(x)

4 Modelización de fenómenos

Al resolver el sistema de forma genérica, no podemos representarlo gráficamente. Es por ello que vamos a fijar valores a algunas variables:

[math] L = 1, \quad \rho = 1, \quad Q = 1, \quad u_0 = 10, \quad u_1 = 1, \quad u_2 = 10 [/math]

De modo que quedan como variables \( K \) y \( \omega(x) \). Además, observamos que:

[math] \alpha = K, \quad f(x) = \omega(x) [/math]

Procedemos a variarlos para modelizar distintas situaciones.

% Ecuación del Calor

clear all; close all;

t0 = 0;         % Instante de tiempo inicial
tf = 1;         % Instante de tiempo final
L = 1;          % Longitud de la barra
a = 1;          % Coeficiente flamígero inicial
alpha = 0.1;    % Coeficiente de difusión
K = 0;          % Temperatura de la fuente de calor externa

u0 = 0;         % Temperatura de la barra en el extremo izquierdo
uL = 0;         % Temperatura de la barra en el extremo derecho

%ui = @(x) a; % Temperatura de la barra en el instante inicial (Constante)
%ui = @(x) a*(L*x - x.^2); % Temperatura de la barra en el instante inicial (Parábola)
ui = @(x) a.*exp(-((x - L/2) / 0.1).^2); % Temperatura de la barra en el instante inicial (Campana de Gauss)

XX = linspace(0,L,500);   % Linspace espacial [0,L]
TT = linspace(t0,tf,500); % Linspace temporal [t0,tf]

v = @(x) K/2*(L*x - x.^2) + (uL - u0)/L*x + u0; % Solución estacionaria (Fuente externa f(x) = K)

% Defino w = u - v

w0 = @(x) ui(x) - v(x); % Condición inicial de w(x,t)

% Cálculo de los coeficientes de Fourier de w(x,t)

n = 10;                   % Nº de términos en la Serie de Fourier
CC = zeros(1,n);          % Vector de coeficientes ck
SS = zeros(1,length(XX)); % Serie de Fourier de w disretizada en espacio

for k = 1:n
        % Aproximación de la integral usando la regla del trapecio
        W = ones(size(XX));  
        W(1) = 1/2; W(end) = 1/2;
        integrand = w0(XX) .* sin(k * pi * XX / L)/sqrt(L);  
        ck = 2 * trapz(XX, integrand .* W);
        CC(k) = ck; % Coeficiente de Fourier c_k
end

w = @(x,t) 0; % Solución del nuevo problema

for k = 1:n % Actualiza w(x)
    w = @(x,t) w(x,t) + CC(k).* exp((-alpha * k^2 * pi^2.* t)/L).* sin(k * pi.* x / L)/sqrt(L);
end

u = @(x,t) w(x,t) + v(x); % Solución del problema original

%% Representación en vídeo 2D %%

MiPeli = VideoWriter('varilla6_2D');  % Crear el vídeo
MiPeli.FrameRate = 30;             % Velocidad de reproducción
open(MiPeli);

for i = 1:length(TT)
    
    figura = figure(1);
    plot(XX,u(XX,TT(i)),color='red',LineWidth=1.5); % Solución en el instante temporal TT(i)
    axis([0 L 0 1])
    xlabel("x"); ylabel("u(x)")
    title("Evolución de la difusión de calor en tiempo")
    grid on
    imagen = getframe(figura);
    writeVideo(MiPeli,imagen); % Insertar imagen en cada frame
    
end

close(MiPeli)


%% Representación en tiempo con vista de planta %%

MiPeli = VideoWriter('Varilla6_VistaPlanta');  % Crear el vídeo
MiPeli.FrameRate = 30;                 % Velocidad de reproducción
open(MiPeli);

for i = 1:length(TT)
    
    figura = figure(1);
    imagesc(u(XX,TT(i)))
    colormap hot
    colorbar
    clim([0 1]);
    axis off
    title("Evolución de la difusión de calor en tiempo")
    imagen = getframe(figura);
    writeVideo(MiPeli,imagen); % Insertar imagen en cada frame
    
end

close(MiPeli)


%% Representación en vídeo 3D %%

MiPeli = VideoWriter('varilla6_3D');  % Crear el vídeo
MiPeli.FrameRate = 30;               % Velocidad de reproducción
open(MiPeli);

for i = 2:length(TT)
    
    figura = figure(1);

    TTi = TT(1,1:i);
    [XXX,TTT] = meshgrid(XX,TTi);
    mesh(XX,TTi,u(XXX,TTT))
    axis([0 L t0 tf])
    xlabel("x"); ylabel("t"); zlabel("u(x,t)")
    title("Representación espacio-temporal de la solución u(x,t)")
    colormap hot
    colorbar

    imagen = getframe(figura);
    writeVideo(MiPeli,imagen); % Insertar imagen en cada frame
    
end

close(MiPeli)

%% Representación tridimensional %%

[XXX,TTT] = meshgrid(XX,TT); % Malla espacio-temporal
mesh(XX,TT,u(XXX,TTT))
axis([0 L t0 tf])
xlabel("x"); ylabel("t"); zlabel("u(x,t)")
title("Representación espacio-temporal de la solución u(x,t)")
colormap hot
colorbar

4.1 Pérdida de calor en la varilla

Suponiendo que no hay fuente de calor externa, es decir f(x) = 0, e imponiendo que la temperatura en ambos extremos sea nula [math] \quad u_1 = 0, \quad u_2 = 0, [/math] observamos la pérdida de calor en la varilla con el paso del tiempo.

Condición inicial de tipo Constante. Alpha = 0.1

Condición inicial de tipo Parabólica. Alpha = 0.1

Condición inicial de tipo Campana de Gauss. Alpha = 0.1

Si suponemos que la temperatura en los extremos de la barra es [math] \quad u_1 = 1, \quad u_2 = 10, [/math] observaremos que se forma un gradiente de temperatura a lo largo de la varilla.

u1 = 1, u2 = 10

4.2 Variación de la constante de difusión [math] \alpha [/math]

Fijados los parámetros anteriores y modificando el coeficiente de difusión [math] \alpha [/math], observamos cómo el calor se difunde más rápido. Físicamente, esto quiere decir que cuanto mayor es la conductividad térmica del material [math] k [/math] y menor es su densidad [math] \rho [/math], más rápido se difunde el calor a lo largo del mismo.

Vista de planta de la varilla. [math] \alpha = 0.1 [/math]

Vista de planta de la varilla. [math] \alpha = 0.5 [/math]

Representación tridimensional de la solución, que avanza en tiempo en el eje Y. [math] \alpha = 0.1 [/math]

Representación tridimensional de la solución, que avanza en tiempo en el eje Y. [math] \alpha = 0.5 [/math]

.

4.3 Aumento de temperatura por una fuente de calor externa

Si la temperatura inicial de la barra es [math] \quad u_0 = 0 [/math] y existe una fuente de calor externa [math] f(x) = K [/math], observamos cómo esta gana calor conforme pasa el tiempo. Con diferentes fuentes comparables, se puede visualizar el Principio de Comparación, que afirma lo siguiente:

Si [math] u,v \in C^{2,1}(Q_T) \cap C(\bar{Q_T}) [/math], donde [math] Q_T [/math] es la región de nuestro dominio hasta tiempo [math] T [/math], son soluciones de la ecuación del calor con fuentes de calor [math] f_1(x), f_2(x) [/math] acotadas en [math] Q_T [/math].

Si

[math] u \geq v [/math] en [math] \partial_pQ_T [/math] y [math] f1(x) \geq f2(x) [/math] en [math] \bar{Q_T} [/math] se tiene que [math] u \geq v [/math] en [math] \bar{Q_T} [/math].

Es decir, si tanto las condiciones fronterizas de ambas soluciones como sus fuentes de calor externas son comparables, se tiene que las soluciones son comparables en cualquier instante de tiempo.

(En este usaría una función f que sea una indicatriz que en un entorno de x=1/3 valga 25 y en un entorno de x=2/3 valga 15, de modo que fisícamente es como si estuviera en contacto en dichos puntos con barras de metal a esas temperaturas (podríamos hacer que una de esas barras estuviera a temperatura negativa a ver que pasa) (podríamos hacer para 2 alphas distintas))

4.4 Modelo 4

(Si que podemos usar una parábola también como función, pues en el caso de la vela la trasmisión de calor se hace por el aire, por lo que si que se transmite calor a toda a barra pero en el centro donde pongas la vela en más medida (debería de ser apuntada))