Placa en forma de Anillo
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Transmisión de calor en placa con forma de anillo. Grupo 11-B |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2013-14 |
| Autores | Javier Chamizo, Jacobo Campos, Miguel García, Javier Mellado, Alfonso Ascanio |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Transmisión de calor en una placa con forma de anillo. Es bastante común modelizar el comportamiento del calor en superficies mediante las ecuaciones diferenciales, y además dada la particularidad de la forma de nuestro caso, utilizando coordenadas polares.
Contenido
1 Introducción
Tenemos la placa que se ve en la imagen superior con las siguientes condiciones. Colocamos en las fronteras interior ρ = 1 y exterior ρ = 6, que sabemos que mantienen una temperatura constante de 0 grados y 10 grados respectivamente como si estuviesen en contacto con un depósito. En nuestro caso vamos a trabajar en coordenadas polares dada la forma de la placa que se presta a ello. Suponemos también que la temperatura u de la placa solo depende de la coordenada radial y el tiempo, es decir u = u(ρ, t). La placa sabemos que satisface la ecuación del calor que se muestra a continuación:
[math] u_t-\Delta u=0[/math]
La placa a su vez dispone de las siguientes condiciones iniciales de temperatura u(p,0):
[math]\left\{\begin{matrix}100·(p-1) → p ∈ [1, 2] \\100 → p ∈ [2, 5]\\90·(6-p)+10 → p ∈ [5, 6]\end{matrix}\right.[/math]
Además tenemos condiciones tipo Dirichlet en los extremos ya que se mantiene constante la temperatura a lo largo del tiempo en estos al tener objetos que mantienen la temperatura constante en contacto con los mismos:
[math]\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.[/math]
Debido a que el ángulo teta y la altura Z se mantienen constantes al ser una figura plana, el Laplaciano en coordenadas cilíndricas es el siguiente:
[math]\Delta u=\frac{1}{p} \frac{\partial }{\partial p} (p \frac{\partial u }{\partial p})+\frac{1 }{p^2} \frac{\partial^2 u }{\partial θ^2}+\frac{\partial^2 u }{\partial z^2}= \frac{\partial^2 u }{\partial p^2}+\frac{1 }{p}\frac{\partial u}{\partial p} [/math]
Al aplicar el Laplaciano, la ecuación del calor que se va a resolver queda de la forma: [math] u_t-u_{pp}-\frac{1}{p} u_p=0 [/math] Quedando así similar al sistema que obtendríamos con una barra con longitud entre 1 y 6.
2 Resolución del Sistema por diferencias finitas
Aquí se va a resolver el problema planteado por el método de diferencias finitas con paso de discretización ∆p = 0.1; y para el tiempo usaremos en cada caso los métodos del trapecio, Euler modificado, explícito e implícito. Todos ellos tomando ∆t = ∆p/4 en t ∈ [0, 10]
2.1 Método del trapecio
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=50;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
%Metodo del trapecio
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%iteraciones desde W^j-->W^j+1
for j=1:length(t)-1
WW=(eye(N-1)+dt/2*(K+L))\((eye(N-1)-dt/2*(K+L))*WW-dt*F);
sol(j+1,:)=[0,WW',10];%guardamos solucion
end
[pp,tt]=meshgrid(p,t);
mesh(pp,tt,sol)
Ahora particularizamos la solución en ρ = 3 y observamos el comportamiento de la temperatura en dichos puntos. La gráfica siguiente muestra una representación 2D temperatura/tiempo
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=50;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
%Metodo del trapecio
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%iteraciones desde W^j-->W^j+1
for j=1:length(t)-1
WW=(eye(N-1)+dt/2*(K+L))\((eye(N-1)-dt/2*(K+L))*WW-dt*F);
sol(j+1,:)=[0,WW',10];%guardamos solucion
end
[pp,tt]=meshgrid(p,t);
J=ones(401,51);
plot3(3*J,tt,sol)
2.2 Método de Euler explícito
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler explicito
WW=WW-dt*((K+L)*WW+F);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
2.3 Método de Euler implícito
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler implicito
WW=(eye(N-1)+dt*(K+L))\(WW-dt*F);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
2.4 Método de Euler modificado
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler explicito
k1=(K+L)*WW+F;
k2=((dt/2)+(K+L))*WW+F;
WW=WW-(dt/2)*(k1+k2);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
3 Evolución de la solución en el tiempo
Se puede deducir a partir de las siguientes gráficas que para grandes tiempos, la temperatura apenas varía dada la solución caracterizada por una función logarítmica decreciente. En este caso el término [math] u_t[/math] se puede despreciar y la temperatura toma un valor estacionario en cada punto de la varilla. Esta ecuación verifica el estado estacionario: [math] u_{pp}+\frac{1}{p} u_p=0 [/math] Resolviendo la ecuación analíticamente observamos que satisface la función de tipo logarítmica [math] u=\frac{10}{log(6)} log(p) [/math], tal como se observa si representamos la solución del problema añadiendo las siguientes lineas a cualquiera de los métodos realizados anteriormente.
plot(p,sol(end,:))
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]';
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold offLa distancia entre la solución estacionaria y la solución u(ρ, t) para t = 0, 1, 2, 10 va aumentando a medida que aumenta el tiempo.
En el caso de dividir la discretización en p por 10 observamos:
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
