Diferencia entre revisiones de «Ecuacion de calor (Grupo 25C)»

Revisión del 07:02 18 may 2015

1 Introducción

Se considera una placa en forma de anillo comprendida entre los radios ρ=1 y ρ=3 y unos objetos colocados en estas fronteras que mantienen la placa a una temperatura constante de 3 y 0 grados respectivamente.

Representación gráfica de la placa

Consideramos también que su temperatura inicial sigue la siguiente función:

centro

Todas estas características llevan al desarrollo de un problema de calor cuyo sistema se va a plantear a continuación.

2 Sistema de ecuaciones en coordenadas polares

En este caso se va a suponer que la temperatura "U" de la placa solamente depende de la coordenada radial y el tiempo, descartándose así la coordenada cíclica theta. Además se sabe que "U" satisface la siguiente ecuación:

centro

Sabiendo esto se puede desarrollar un sistema de ecuaciones en coordenadas polares que representa el problema planteado.

Al no tener dependencia de la variable “θ” el problema que nos queda es de contorno en “ρ”:

[math] \Delta U = U_{ρρ} +\frac{1 }{ρ} U_ρ; (1) [/math]

Aplicando la ecuación que nuestro problema tiene que satisfacer:

[math] U_t - \Delta U = 0; (2)[/math]

Sustituyendo (1) en (2):

[math] U_t - U_{ρρ} -\frac{1}{ρ} U_ρ = 0 [/math]

Añadiendo las condiciones de contorno e inicial nos queda el siguiente planteamiento:

[math] \left\{\begin{matrix} U_t - U_{ρρ} -\frac{1 }{ρ} U_ρ = 0 \\ U(1,t)=3 ; U(3,t)=0 \\U(ρ,0)=e^{-8(ρ-2)^2}+3-ρ \\ ρ∈(1,3);t\gt0 \end{matrix}\right.[/math]

3 Resolución del problema

Se resuelve por el método de diferencias finitas discretizando la variable espacial "ρ" y la variable temporal "t". En este caso tenemo dos matrices K1 y K2 que acompañan el vector "U" ya que el planteamiento lo exige por las apariciones de las derivadas primera y segunda de U con respecto de "ρ".

3.1 Método del trapecio

En este caso se pide realizar el metodo con los siguientes pasos:

clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);Ua=3;Ub=0;
F(1)=F(1)+Ua*((1/h^2)-(1/(2*h*rr(1))));
F(end)=F(end)+Ub*((1/h^2)+(1/(2*h*rr(end))));
%Discretizacion temporal
t0=0;tM=10;
k=h/4;
t=t0:k:tM;
M=length(t)-1;
%Definicion del valor inicial
U0=exp(-8.*((rr-2).^2))+3-rr;
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+(k/2)*(K1+K2))\((eye(size(K1))-(k/2)*(K1+K2))*sol(:,i)+k*F');
end
%Grafico
UA=Ua*ones(1,length(t));
UB=Ub*ones(1,length(t));
sol=[UA;sol;UB];
[Mr,Mt]=meshgrid(r,t);
s=surf(Mr,Mt,sol');
set(s,'edgecolor','flat');
colormap summer

Gráfico del problema mediante el método del trapecio

3.2 Temperatura para ρ=2

Para comprobar como cambia la temperatura en función del tiempo para ρ=2 hay que implementar el siguiente codigo en Matlab, a continuación del anterior:

d=(2-a)/h+1;
plot(t,sol(d,:));

Como resultado nos sale el siguiente gráfico donde se observa que para un radio 2, la temperatura baja desde los dos grados hasta los 1.1 para luego mantenerse constante:

Gráfico para ρ=2

3.3 Método de Euler explícito

Para obtener tanto este método como el de Euler implícito y el de Heun, lo unico que hay que modificar es el apartado "%Solucion" del código del método del trapecio.

%Solucion
for i=1:M
    sol(:,i+1)=sol(:,i)+k*(-(K1+K2)*sol(:,i)+F');
end

En este caso se han tenido problemas con la obtención total de los datos y del gráfico para un tiempo de 10 segundos ya que el método es muy impreciso y da valores anormalmente altos. Para un t=1 se ha podido obtener el siguiente gráfico:

Gráfico con el método de Euler explítcito

3.4 Método de Euler implícito

%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end

La gráfica obtenida es muy parecida a la del método del trapecio:

Gráfico meiante el método de Euler implícito

3.5 Método de Heun

El código empleado ha sido el siguiente:

%Solucion
for i=1:M
    Q1=(-(K1+K2))*sol(:,1)+F';
    Q2=(-(K1+K2))*(sol(:,1)+k*Q1)+F';
    sol(:,i+1)=sol(:,1)+(k/2)*(Q1+Q2);
end

Gráfico mediante el método de Heun

4 Planteamiento del problema para tiempos grandes

En este caso es el mismo planteamiento que el anterior pero sin la variable dependiente del tiempo. El problema se reduce a un problema de contorno que se resuelve mediante el siguiente código:

clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);
Ua=3;Ub=0;
F(1)=F(1)+Ua*((1/h^2)-(1/(2*h*rr(1))));
F(end)=F(end)+Ub*((1/h^2)+(1/(2*h*rr(end))));
%Discretizacion temporal
Y=(K1+K2)\F';
Y=[Ua;Y;Ub];
plot(r,Y);

La gráfica resultante nos muestra a qué debería parecerse la función de temperatura para tiempos grandes:

centro

Volviendo a plantear el sistema de antes con Euler implícito y añadiendo los valores temporales para los cuales queremos estudiar el error de los resultados

clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);
Ua=3;Ub=0;
F(1)=F(1)+Ua*((1/h^2)-(1/(2*h*rr(1))));
F(end)=F(end)+Ub*((1/h^2)+(1/(2*h*rr(end))));
%Discretizacion temporal
t0=0;tM=10;
k=h/4;
t=t0:k:tM;
M=length(t)-1;
%Definicion del valor inicial
U0=exp(-8.*((rr-2).^2))+3-rr;
Y=(K1+K2)\F';
Y=[Ua;Y;Ub];
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end
%Grafico
UA=Ua*ones(1,length(t));
UB=Ub*ones(1,length(t));
sol=[UA;sol;UB];
Y1=sol(:,1);
Y2=sol(:,41);
Y3=sol(:,81);
Y4=sol(:,401);
figure(1)
hold on
plot(r,Y,'r');
plot(r,Y1,'g');
plot(r,Y2,'b');
plot(r,Y3,'m');
plot(r,Y4,'k');
legend('Estacionaria','Para t=0','Para t=1','Parat=2','Para t=10');
hold off
error1=abs(Y(2)-Y1(2));
error2=abs(Y(2)-Y2(2));
error3=abs(Y(2)-Y3(2));
error4=abs(Y(2)-Y4(2));
errores=[error1;error2;error3;error4];

El gráfico resultante nos muestra que a medida que avanza el tiempo el error es cada vez mas pequeño

Gráfico de errores

Tomando el punto central ρ=2, estos son los valores obtenidos de los errores

5 Planteamiento del problema con una pieza aislante en ρ=1

En este caso el problema cambia ligeramente y obliga a que la condicion de contorno en el interior del disco sea:

[math] U_ρ(1,t)=0 [/math]

Siguiendo el procedimiento y aproximando la derivada parcial:

[math] U_ρ = \frac{U_{n+1}(t)-U_{n-1}(t)}{2h} = 0 [/math]

Donde para el radio interior, es decir n=0 la ecuación queda en:

[math] U_1(t)=0 [/math]

Se ha utilizado el siguiente código para resolver este problema:

clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K1(1,1)=-1;K1(2,1)=0;
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
K2(2,1)=0;K2(1,1)=1;
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);Ub=0;
%Discretizacion temporal
t0=0;tM=10;
k=h/4;
t=t0:k:tM;
M=length(t)-1;
%Definicion del valor inicial
U0=exp(-8.*((rr-2).^2))+3-rr;
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end
A=zeros(1,length(t));
%Grafico
UB=Ub*ones(1,length(t));
sol=[A;sol;UB];
sol(1,:)=sol(2,:);
sol(2,:)=0;
[Mr,Mt]=meshgrid(r,t);
s=surf(Mr,Mt,sol');
set(s,'edgecolor','flat');
colormap summer

En el grafico resultante se puede observar como la temperatura del radio interior sube a medida que avanza el tiempo debido al aislamiento térmico

Gráfico para ρ=1 aislado

En cuanto al valor estacionario se puede deducir analíticamente que es U=0.

6 Disco de radio 3

En este apartado consideramos un disco completo de radio 3 que cumple la misma ecuación que el disco anterior con la diferencia de que en este caso la condicion inicial es: [math] u(3,t)=0 [/math]

6.1 Problema de autovalores

La condicion de que la funcion sea diferenciable en el punto ρ=0 nos lleva a la conclusion de que la derivada en ese mismo punto es igual a 0:

[math] u_ρ(0,t)=0 [/math]

A partir de este resultado y con la ecuación inicial, tenemos el siguiente sistema:

[math] \left\{\begin{matrix} U_t - U_{ρρ} -\frac{1 }{ρ} U_ρ = 0 \\ U_ρ(0,t)=0 ; U(3,t)=0 \\ t\gt0 \end{matrix}\right.[/math]

Aplicando el método de separación de variables donde U(ρ,t)=ρT llegamos al desarrollo de dos problemas, uno de contorno en ρ y uno de valor inicial en T. El problema de contorno en ρ queda de la siguiente manera:

[math] \left\{\begin{matrix} ρ'' + \frac{1}{ρ}ρ'+λρ=0 \\ ρ'(0)=0;ρ(3)=0 \end{matrix}\right.[/math]

6.2 Bessel

Tomando la ecuacion anterior y multiplicándola por [math]ρ^2[/math] da como resultado una nueva ecuación igual que la que se presenta en el ejemplo del enunciado:

[math] \left\{\begin{matrix} {ρ^2}ρ'' + {ρ}ρ'+{ρ^2}λρ=0 \\ ρ'(0)=0;ρ(3)=0 \end{matrix}\right. [/math]

Las autofunciones [math] φ_k = J_0 (ρ{λ}^{1/2}) [/math] son soluciones de este problema y se puede demostrar derivando y sustituyendo en el problema que queremos resolver. Los ceros de la funcion de Bessel son los valores para los cuales la funcion es nula. Por lo tanto, en este caso tienen que ser los valores [math] ρ{λ}^{1/2} [/math] para los cuales las autofunciones [math] φ_k = J_0 (ρ{λ}^{1/2}) [/math] son cero. Tomamos [math] C_0 [/math] a los ceros de estas autofunciones:

[math] C_0 = ρ{λ}^{1/2} [/math]

Ahora, aplicamos la condicion de contorno [math] ρ(3)=0 [/math] , de lo que resulta:

[math]φ_k = J_0 (3{λ}^{1/2}) = 0 [/math] [math]C_0 = 3{λ}^{1/2}[/math]