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 | |
Contenido
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 ut-Δu=0 , siendo Δu=0 el laplaciano de u, 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');
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.
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
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
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
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 ut, de esta manera la temperatura toma un valor estacionario en cada punto de la placa. En este caso la ecuación quedaría así: -utt-1/ρut=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
%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 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.
Para apreciar la diferencia entre la solución estacionaria y la solución u(ρ,t) para t=0,1,2,10 se usa el programa anterior estudiando los errores entre ellos.
%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,'k');
plot(r,Y1,'b');
plot(r,Y2,'r');
plot(r,Y3,'g');
plot(r,Y4,'p');
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];
