Ecuación del calor. (Grupo A8)

Trabajo realizado por estudiantes
Título	Ecuación del calor (Grupo A8)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Valentina Salazar; Antonio Carrero; José Francisco Aguilera
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción y modelización del problema.

En este artículo trataremos el comportamiento térmico de una varilla sometida a ciertas condiciones térmicas y físicas mediante diferentes métodos numéricos, daremos interpretación física a los resultados arrojados por estos métodos. Basaremos el estudio analítico y numérico necesario en la ecuación del calor propuesta por Jean-Baptiste Joseph Fourier en 1822.

Para estudiar el problema consideraremos una varilla delgada, de sección constante y de un material homogéneo, de longitud [math]L=4[/math]. La situaremos en el intervalo [math]x\in{(0,4)}[/math] de la recta real. Para llevar acabo esta modelización tendremos que asumir unas ciertas hipótesis y simplificaciones:

• Para el primer caso a plantear la varilla estará infitamente aislada del entorno y por tanto no habrá flujo de calor sobre la superficie lateral de la misma.
• Al estar considerando el caso unidimensional de la ecuación del calor, asumiremos que al ser una varilla delgada la temperatura a lo largo de una sección ortogonal al eje [math]x[/math] se mantiene constante en toda la sección.
• Consideraremos que el calor específico del material, [math]c[/math], es constante y no depende de la temperatura, y por tanto la difusividad térmica,[math]\alpha=\frac{k}{c\rho}[/math] también lo será.

Si damos un valor a las constantes [math]c[/math], [math]k[/math] y [math]\rho[/math], tal que [math]\alpha=2[/math]; y no hay fuentes ni sumideros de calor dentro de la varilla, la ecuación en derivadas parciales que tendrá que cumplir la distribución de temperaturas a lo largo de la misma es::

[math]u_t(x,t)-2u_{xx}(x,t)=0[/math]

Donde:

• [math]u(x,t)[/math] nos define como evoluciona la temperatura a lo largo del eje [math]x[/math].
• [math]u_t(x,t)[/math] es la derivada de la temperatura con respecto del tiempo.
• [math]u_{xx}(x,t)[/math] es la segunda derivada de la temperatura con respecto de [math]x[/math].

Junto con está ecuación en derivadas parciales consideraremos una serie de condiciones: iniciales y de contorno. Como condición inicial para está modelización propondremos::

[math]u(x,0)=e^{-6(x-2)^2}+5-\frac{5}{4}x[/math]

Donde [math]u(x,0)[/math] no es otra cosa que la distribución de temperaturas inicial.

A lo largo del artículo variaremos las condiciones de contorno para dar diferentes ejemplos, pero en el bloque de casos prácticos trabajaremos con las siguientes::

[math]u(0,t)=5[/math]: [math]u(4,t)=0[/math]

Éstas indican la temperatura (en este caso, otras pueden ser los flujos de calor o combinaciones de ambos) en cada uno de los extremos de la varilla.

2 Casos prácticos

En los apartados que siguen estudiaremos a través del cálculo numérico el caso anteriormente expuesto e interpretaremos los resultados obtenidos mediante dicho análisis. El lenguaje usado en el que se implementarán estos modelos será código Matlab u Octave.

2.1 Método de diferencias finitas o método de lineas

Trataremos de reducir el sistema de ecuaciones a una sola variable en una discretización de [math]n[/math] valores (siendo [math]h[/math] el tamaño de los subintervalos) haciendo la aproximación::

[math]u_{xx}(x,t)=\frac{-U^{N-1}(t)+2U^N(t)-U^{N+1}(t)}{h^2}[/math]

Así podemos aplicar uno de los diferentes métodos iterativos en la variable [math]t[/math] que se muestran en los subapartados que siguen.

2.1.1 Método del trapecio

Se trata de un método implícito, es decir que hay que despejar en la ecuación matricial que nos da la [math]U(t)[/math]. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
    sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)
 
 
 
[xx,tt]=meshgrid(x,t);
 
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')

Obteniendo los siguientes gráficos:

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.

Como podemos ver la temperatura apenas tarda en estabilizarse en el tiempo a lo largo de los puntos de la varilla, esto nos sugiere que una buena aproximación será la solución estacionaria de la ecuación.

Sucesivas temperaturas del punto medio frente al tiempo.

La temperatura del punto medio descenderá hasta volverse casi asintótica hacia un valor aproximado de 2.625.

2.1.2 Método de Euler explícito

A diferencia del método del trapecio, como su nombre indica, este método es explicito, lo que es una ventaja a la hora de implementar ya que nos ahorramos el despejar; pero sin embargo se trata de un método bastante inestable, y para que arroje soluciones lógicas hemos de disminuir el tamaño del paso temporal frente al espacial del orden de [math]dt=\frac{h^2}{4}[/math], siendo [math]h[/math], [math]dt[/math] dichos pasos espacial y temporal respectivamente. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
    sol(:,n+1)=M*sol(:,n)+dt*F;
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion
 
 
 
[xx,tt]=meshgrid(x,t);
 
figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')

Obteniendo los gráficos: