Revisión actual del 01:03 20 mar 2025

1 Introducción y enfoque

La modelización matemática se emplea en 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.

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.
• 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.

2 Planteamiento del sistema de EDP

Imaginemos que tenemos una varilla metálica con densidad constante \( \rho \), de longitud \( L \), cuya temperatura inicial es \( u_0(x) \). Además, los extremos de la varilla mantienen temperaturas fijas \( u_1 \) y \( u_2 \).

También consideramos una función \( \omega \), que representa la producción de energía externa, la cual supondremos constante en tiempo y espacio.

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 [/math]

donde \( k \) es la conductividad térmica de la varilla y \( Q \) es el calor específico, que relaciona la temperatura con la energía.

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

Con esto, si aplicamos nuestras condiciones iniciales y fronterizas, el sistema queda:

[math] \begin{cases} u_t - \alpha u_{xx} = f, & x \in [0,L], \quad t \gt 0 \\ u(x,0) = u_0(x), & 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

Suponiendo que [math] f = K [/math] constante, la solución estacionaria es [math] v(x) = \dfrac{K}{2\alpha}\left(Lx -x^2 \right) + \dfrac{(u_2 - u_1)}{L} x + u_1 [/math]. Definiendo [math] w(x,t) = u(x,t) - v(x) [/math] y aplicando separación de variables, obtenemos

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

y [math] u(x,t) = \sum_{n=1}^{\infty} \dfrac{c_n}{\sqrt{L}} e^{- \frac{n^2\pi^2}{\alpha L} t} \sin{\left( \dfrac{n\pi x}{L} \right)} + v(x) [/math] es la solución del problema original.

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_1 = 0, \quad u_2 = 0 [/math]

Y variamos los demás parámetros para modelizar distintas situaciones con el siguiente código en MatLab

% Ecuación del Calor

clear all; close all;

t0 = 0;         % Instante de tiempo inicial
tf = 0.5;         % 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*alpha)*(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^2).* 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 [math] f(x) = 0 [/math] y con temperatura nula en ambos extremos [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. [math] \alpha = 0.1 [/math]

Condición inicial de tipo parabólica. [math] \alpha = 0.1 [/math]

Condición inicial de tipo campana de Gauss. [math] \alpha = 0.1 [/math]

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

u1 = 1, u2 = 10

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

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 varilla es [math] u_0 = 0 [/math] y existe una fuente de calor externa [math] f(x) = K \gt0 [/math], observamos cómo esta gana calor conforme pasa el tiempo. Con diferentes fuentes comparables se cumple el principio de comparación, que afirma:

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] \begin{cases} u(x,t) \geq v(x,t) \text{ en } \partial_pQ_T \\ f_1(x) \geq f_2(x) \text{ en } \bar{Q_T} \end{cases} [/math]

Se tiene que [math] u \geq v [/math] en [math] \bar{Q_T} [/math].

Es decir, si las condiciones fronterizas de ambas soluciones y sus fuentes de calor externas son comparables, las soluciones son comparables.

Vista de planta de la varilla. [math] K = 1 [/math]

Vista de planta de la varilla. [math] K = 3 [/math]

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

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

Al aumentar la longitud de la varilla [math] L [/math] observamos cómo el máximo de la temperatura que alcanza la barra es ligeramente mayor y que el intervalo donde se alcanza es mayor. Mientras que si disminuimos [math] L [/math], se observa cómo al mantenerse los extremos temperatura 0, y ahora estar todos los puntos intermedios más próximos a estos extremos, se llega a reducir el máximo que se alcanza.

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

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