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.
Aparece un término más de lo habitual, la derivada primera en ρ, por lo que habrá que hacer modificaciones con respecto a la resolución típica de estos problemas.
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 líneas con la siguiente discretización: Δρ=0.1, Δt=Δρ/4 en t∈[0,10]. Las condiciones de frontera son ambas Dirichlet. 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. Además se puede comprobar que en ρ=2 la gráfica da una temperatura inicial de 2 grados, lo cual también se obtiene sustituyendo ρ=2 en la fórmula de la temperatura 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 bruscamente hasta 1.1 grados y mantenerse asintóticamente en ese valor.
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 líneas. 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. Por lo tanto, es díficil hacer valoraciones sobre el comportamiento de la temperatura, ya que el resultado obtenido no es en absoluto fiable.
4.2 Método de Euler Explícito
Se usa de nuevo el método de líneas 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.
Cabe destacar como cambian las gráficas de los métodos explícitos cambiando el tamaño de paso a uno mayor, en este caso h=1.
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í: -uρρ-1/ρuρ=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. La gráfica resultante en 2D es una función de la temperatura que sólo depende de la coordenada axial.
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];
Para el siguiente gráfico se ha cambiado el tamaño de paso de ρ a 0.01.
Se aprecia que según aumenta el tiempo se aproxima más la función a la solución estacionaria, como se había comprobado antes que para tiempos grandes la temperatura no varía prácticamente nada.
Estudiando el punto ρ=2 se obtienen las siguientes distancias: con Δρ=0.1 con Δρ=0.01
en t=0, 0.8383 en t=0, 0.9824
en t=1, 0.0080 en t=1, 0.0011
en t=2 7.79e-04 en t=2, 0.0033
en t=10, 6.52e-12 en t=10, 7.80e-04
Esto nos demuestra que las distancias con la función estacionaria, al disminuir el tamaño de paso, disminuyen para tiempos grandes ya que en esos tiempos la aproximación de la solución es mejor a la función u(ρ,t) y a la función estacionaria.
6 Aislante en la frontera interior
En este caso, en lugar de mantener la temperatura a 3 grados como se hacía antes, se coloca una pieza aislante en la frontera interior de la pieza (ρ=1). Con esto se consigue que el flujo con el exterior por ese borde sea nulo. Esto supone un cambio en el problema, ya que ahora una de las condiciones frontera es tipo Neumann. El nuevo problema a resolver es el siguiente:
Por lo que ahora se incluyen las modificaciones que supone la condición Neumann y se resuelve numéricamente con el método del trapecio
% 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) 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
% Ejerciocio tipo Condiciones Dirichlet
rr=r(1: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))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1));
K(1,2)=-2*q/h^2;
L=(1/(2*h))*(-diag(ones(1,N-1),1)+diag(ones(1,N-1),-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,M+1);
sol(:,1)=U0;
% Aplicamos el metodo de Trapecio
for i=1:M
G1=g(rr,t(i));
G1(1)=G1(1)-2*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)-2*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
Ub=cb(t);
sol=[sol;Ub]; % las añadimos
% gráfico
[Mt,Mr]=meshgrid(t,r);
figure (1)
surf(Mr,Mt,sol)
xlabel('rho');
ylabel('t');
figure (2)
plot(t,sol(11,:))
Se ve claramente en la figura como la temperatura estacionaria es u=0. Se alcanza dicho valor a partir de los 8 segundos.
