Ecuación del calor en placa con forma de anillo

Trabajo realizado por estudiantes
Título	Ecuación del Calor en una placa en forma de anillo (Grupo 12)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2016-17
Autores	Nicolás de la Fuente Ceñal Fabio Torres Salas Álvaro Solís González
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Descripción del problema

Se considera una placa plana en forma de anillo comprendida entre los radios ρ=1, ρ=3.Inicialmente la temperatura viene dada por:

u(ρ,0)=e-8(ρ-2)^2-ρ+3.

Se colocan en las fronteras interior y exterior objetos que mantienen una temperatura constante de 3 y 0 grados respectivamente.

Se supone que la temperatura de la placa sólo depende de la coordenada radial y del tiempo, u=u(ρ,t), y que satisface la ecuación del calor u_t-Δu=0 , siendo Δu=0 el laplaciano de u, en polares.

Desarrollo del laplaciano en polares

2 Resolución del problema con el método del Trapecio

Al tratarse de un método implícito, hay que despejar en el esquema numérico, obteniendo el bucle que se mostrará en el programa. Para resolver problema se emplea el método de diferencias finitas con la siguiente discretización: Δρ=0.1, Δt=Δρ/4 en t∈[0,10] Hay que destacar que a diferencia de la programación habitual, se incluye el término de u_ρ. Esto implica que se tenga que añadir la matriz L al desarrollo numérico. U'=-(K+L)*U+G

% Datos del problema

t0=0;
tM=10; %instante final
q=input(' Introduce difusividad');  %1  % difusividad= conductividad termica/densidad y calor especifico

g=@(r,t) 0*r;
ca=@(t) 3+0*t;
cb=@(t) 0*t;
u0=@(r) exp(-8*(r-2).^2)+3-r;

% discretizacion espacial
h=0.1;%Anchura de paso en el espacio
N=round((3-1)/h); %subintervalos en x
r=linspace(1,3,N+1);
r=r'; % El vector r tiene N+ 1 elementos

% Ejercicio tipo Condiciones Dirichlet

rr=r(2:end-1); % tomamos los elementos del segundo al penultimo (N-1 ecuaciones)(N componentes)

%Matriz K de N-1 x N-1)
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));

L=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
A=diag(1./rr);
L=A*L;

%Valor inicial
 U0=u0(rr);
 
% Vector de terminos independientes 
% depende del tiempo ---> queda para despues

k=h/4; 
M=round((tM-t0)/k); % numero de subintervalos en t
t=linspace(t0,tM,M+1); % vector con los tiempos

% sol matriz que tendra la solucion para cada punto e intante
sol=zeros(N-1,M+1);
sol(:,1)=U0;

% Aplicamos el metodo de Trapecio
for i=1:M
    G1=g(rr,t(i));
    G1(1)=G1(1)+q*ca(t(i))/h^2-(ca(t(i)))/(2*h*rr(1));
    G1(end)=G1(end)+q*cb(t(i))/h^2+(cb(t(i)))/(2*h*rr(end));
     G2=g(rr,t(i+1));
    G2(1)=G2(1)+q*ca(t(i+1))/h^2-(ca(t(i)))/(2*h*rr(1));
    G2(end)=G2(end)+q*cb(t(i+1))/h^2+(cb(t(i)))/(2*h*rr(end));
    sol(:,i+1)=(eye(size(K))+(k/2)*(K+L))\(sol(:,i)+k/2*(-(K+L)*sol(:,i)+G1+G2));
end

% tenemos la solucion completa SALVO las condiciones DIRICHLET
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];  % las añadimos

% gráfico

[Mt,Mr]=meshgrid(t,r);
     
surf(Mr,Mt,sol) 

xlabel('x');
ylabel('t');

Gráfica con el método del trapecio

Se comprueba que la gráfica tiene valor 3 en ρ=1 y valor 0 en ρ=3, que son los valores en los extremos para la tempertura inicial.

3 Temperatura con ρ=2

Para observar el comportamiento de la temperatura en ρ=2 se incluye el siguiente comando al final del programa previo

plot(t,sol(11,:))

Con ello se obtiene una gráfica que muestra la temperatura en función del tiempo en el radio=2. Se observa como la temperatura inicial es de 2 grados, como es esperable si entramos con ρ=2 en la función de la temperatura inicial; para luego descender rápidamente con el tiempo.

Comportamiento de la temperatura en los puntos de rho=2 en función del tiempo.

4 Resolución del problema con Heun, Euler y Euler implícito

4.1 Método de Heun

Se vuelve a utilizar el método de diferencias finitas. En este caso se sustituye el bucle por esquema numérico de Heun que es el siguiente

for i=1:M
    G1=g(rr,t(i));
    G1(1)=G1(1)+q*ca(t(i))/h^2-ca(t(i))/(2*h*rr(1));
    G1(end)=G1(end)+q*cb(t(i))/h^2+cb(t(i))/(2*h*rr(end));
     G2=g(rr,t(i+1));
    G2(1)=G2(1)+q*ca(t(i+1))/h^2-ca(t(i))/(2*h*rr(1));
    G2(end)=G2(end)+q*cb(t(i+1))/h^2+cb(t(i))/(2*h*rr(end));
    K1=-K*sol(:,i)+G1;
    K2=-K*(sol(:,i)+K1*k)+G2;
    sol(:,i+1)=sol(:,i)+(k/2)*(K1+K2);
end

Gráfica con el método de Heun

El método de Heun, al contrario que el trapecio, es explícito. Por ello es inestable y se aprecia en la gráfica que tiene picos muy elevados.

4.2 Método de Euler Explícito

Se usa de nuevo el método de diferencias finitas con la introducción del esquema numérico correspondiente al Método de Euler Explícito:

for i=1:M
    Gi=g(rr,t(i));
    Gi(1)=Gi(1)+q*ca(t(i))/h^2-(ca(t(i)))/(2*h*rr(1));
    Gi(end)=Gi(end)+q*cb(t(i))/h^2+(cb(t(i)))/(2*h*rr(end));
    sol(:,i+1)=sol(:,i)+k*(-(K+L)*sol(:,i)+Gi);
end

Gráfca para el método de Euler explícito

Con este método vuelve a ocurrir lo mismo, al ser explícito es inestable y por ello da esos valores tan irregulares.

4.3 Método de Euler implícito

Al tratarse de un método implícito hay que despejar el esquema numérico como se hizo con el del trapecio.

for i=1:M
    Gi=g(rr,t(i));
    Gi(1)=Gi(1)+q*ca(t(i))/h^2-(ca(t(i)))/(2*h*rr(1));
    Gi(end)=Gi(end)+q*cb(t(i))/h^2+(cb(t(i)))/(2*h*rr(end));
    sol(:,i+1)=(eye(size(K))+k*(K+L))\(sol(:,i)+k*Gi);
end

Gráfica euler implícito

Al ser un método implícito, no es inestable y por eso salen resultados más razonables. De hecho es la misma gráfica que en trapecio.

5 Evolución con el tiempo

Se aprecia en las gráficas de Trapecio y Euler implícito que la temperatura apenas cambia para tiempos grandes. Por ello se puede despreciar u_t, de esta manera la temperatura toma un valor estacionario en cada punto de la placa. En este caso la ecuación quedaría así: -u_tt-1/ρu_t=0 Lo cual es una ecuación diferencial de segundo orden que ahora sólo depende de la dirección radial. A ella se añaden las condiciones de contorno u(1)=3; u(3)=0

centro

%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)

Gráfica de la temperatura en función de la coordenada axial

La solución u(ρ,t) debería parecerse a la función u(ρ,0)=e^{-8(ρ-2)^2}-ρ+3. Ya que esta función no depende del tiempo y muestra la temperatura inicial de la placa, la cual se ha probado que no varía para tiempos grandes.