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.

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

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.