Revisión del 22:58 14 may 2015

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:

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

Al igual que en el apartado anterior vemos que se estabiliza muy rápido. Debido al tener que disminuir bastante el tamaño del paso temporal hemos tenido que cambiar la función de Octave para representar este gráfico.

Sucesivas temperaturas del punto medio frente al tiempo.

Vemos que en esta implementación los cambios de temperatura en el punto medio son más continuos, pero también acaba siendo asintótica en un valor aproximado de 2.625.

2.1.3 Método de Euler implícito

Es más fácil hacerlo converger que su homónimo explícito, pero tiene la desventaja de tener que despejar en la ecuación para implementarlo. A continuación se muestra el código que lo implementa en Matlab u Octave, dándonos 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*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)
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.

Volvemos a ver que la estabilización de temperaturas en la varilla es bastante rápida.

Sucesivas temperaturas del punto medio frente al tiempo.

Una vez más observamos un descenso asintótico hacia el valor de temperatura 2.625.

2.1.4 Método de Runge-Kutta

Se trata de un método explícito, por tanto necesitaremos aumentar el tamaño del paso una vez más hasta los mismos valores que para Euler explícito para asegurar su convergencia. 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
   k1=-K*sol(:,n)+F;
   k2=-K*(sol(:,n)+k1*dt/2)+F;
   k3=-K*(sol(:,n)+k2*dt/2)+F;
   k4=-K*(sol(:,n)+k3*dt)+F;
   sol(:,n+1)=sol(:,n)+(dt/6)*(k1+2*k2+2*k3+k4);
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 las siguientes gráficas:

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

Al igual que en todos los métodos anteriores se produce la ya esperada estabilización de temperaturas en la varilla de manera bastante acelerada.

Sucesivas temperaturas del punto medio frente al tiempo.

De nuevo aparece un descenso asintótico hasta el valor 2.625.

2.2 Método de Fourier

Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Los autovalores y autovectores obtenidos para este problema y con condiciones de contorno son::

[math]\mu_k=\frac{k^2\pi^2}{16}[/math]:

[math]\varphi_k(x)=sen(\frac{k\pi x}{4})[/math]:

[math]k=1,2,3,...[/math]

Dados estos valores y la función de homogenización para las condiciones de contorno::

[math]v(x,t)=5-\frac{5}{4}x[/math]

Mediante el siguiente código de Matlab implementamos el problema numéricamente, obteniendo así la gráfica en 3D de la solución para [math]k=20[/math] (por ser la más representativa y exacta),y la comparaciónde las gráficas en el plano de [math]t=0,5[/math] para [math]k=1,3,5,10 y 20[/math].

clear all
N=100;
a=0; b=4;
h=(b-a)/N; x=0:h:4; % space discretization
ht=h; t=0:ht:5; % time discretization
[xx,tt]=meshgrid(x,t); % time discretization
f=(exp((-6)*(x-2).^2)); % Initial data
aprox=0; % Compute the approximation
q=1;
for k=1:q
  autovalor=((k*pi)^2)/16; % eigenvalue
  autofuncion=sin(k*pi*x/4); % eigenfunction
  c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
  aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
  end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:)) % Draw the approximation
hold on
q=3;
for k=1:q
  autovalor=((k*pi)^2)/16; % eigenvalue
  autofuncion=sin(k*pi*x/4); % eigenfunction
  c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
  aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
  end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'r') % Draw the approximation
hold on
q=5;
for k=1:q
  autovalor=((k*pi)^2)/16; % eigenvalue
  autofuncion=sin(k*pi*x/4); % eigenfunction
  c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
  aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
  end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'m') % Draw the approximation
hold on
q=10;
for k=1:q
  autovalor=((k*pi)^2)/16; % eigenvalue
  autofuncion=sin(k*pi*x/4); % eigenfunction
  c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
  aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
  end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'g') % Draw the approximation
hold on
q=20;
for k=1:q
  autovalor=((k*pi)^2)/16; % eigenvalue
  autofuncion=sin(k*pi*x/4); % eigenfunction
  c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
  aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
  end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'y') % Draw the approximation
hold on
plot(x,5-5/4*x,'k')
figure(2)
mesh(xx,tt,aprox2)
hold on

Obteniendo así:

Gráfico que muestra de la temperatura de la varilla frente a al tiempo a lo largo de sus diferentes posiciones.

Se puede ver como a diferencia del método de lineas el error viene inducido por la aproximación mediante series de Fourier de la condición inicial. Por lo demás evoluciona de manera similar a los métodos aplicados anteriormente.

Gráficas de las temperaturas a lo largo de la varilla con aproximaciones de diferentes índices del sumatorio en la solución. k=1(azul). k=3(rojo). k=5(magenta). k=10(verde). k=20(amarillo). solución estacionaria(negro)

En la gráfica anterior se puede apreciar el nivel de aproximación con relación al número de coeficientes que utilizamos de la serie. Vemos como cuantos menos coeficientes usamos más nos acercamos incluso en tiempo cortos a la solución estacionaria, lo que nos indica que desciende el nivel de precisión de estas aproximaciones.

3 Cambio de condiciones de contorno

Ahora manteniendo la ecuación que nos modelizaba el comportamiento térmico de la varilla, cambiamos las condiciones de contorno en los extremos. Las nuevas condiciones propuestas se pueden asimilar como el aislamiento térmico total del extremo izquierdo, lo que anula el flujo de calor en ese lado. Estas son:

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

Para implementarlo en sus nuevas condiciones usaremos el método de lineas con el método del trapecio. El código que implementa esta modelización en lenguaje Matlab y nos proporciona un gráfico de la temperatura frente a las variables espacial y temporal es:

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=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,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);
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
UB=zeros(length(t),1);
sol=[sol;UB'];
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')

4 Interpretaciones estacionarias

5 Intercambio de calor con el entorno. Varilla no aislada longitudinalmente

5.1 Cambio de condiciones de contorno