Diferencia entre revisiones de «Placa en forma de Anillo»

De MateWiki
Saltar a: navegación, buscar
 
(No se muestran 35 ediciones intermedias de 2 usuarios)
Línea 7: Línea 7:
 
[[Archivo:placadescripcion.png|marco|derecha|Placa en anillo con p ∈ [1, 6]]]
 
[[Archivo:placadescripcion.png|marco|derecha|Placa en anillo con p ∈ [1, 6]]]
  
−
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. La placa a su vez dispone de las siguientes condiciones iniciales de temperatura u(p,0):
+
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>
 
<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>
Línea 15: Línea 19:
 
<math>\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.</math>
 
<math>\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.</math>
  
−
Por otro lado la ecuación del calor es:
+
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:
  
−
Sin embargo al tomar el Laplaciano esta queda transformada en:
+
<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>
−
 
+
−
Ya que el ángulo teta y la altura Z se mantienen constantes 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.
 
Quedando así similar al sistema que obtendríamos con una barra con longitud entre 1 y 6.
  
Línea 29: Línea 31:
 
===Método del trapecio===
 
===Método del trapecio===
 
{{matlab|codigo=
 
{{matlab|codigo=
 +
%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)
 
}}
 
}}
 +
[[Archivo:trapeciobueno.png|marco|centro|Distribución de la temperatura en el disco]]
 +
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
 +
{{matlab|codigo=
 +
%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)
 +
}}
 +
 +
[[Archivo:p3.png|marco|centro|Distribución de la temperatura para p=3]]
 +
 
===Método de Euler explícito===
 
===Método de Euler explícito===
 
{{matlab|codigo=
 
{{matlab|codigo=
 +
%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)
 
}}
 
}}
 +
 +
[[Archivo:explicitobueno.png|marco|centro|Distribución de la temperatura en el disco]]
 +
 
===Método de Euler implícito===
 
===Método de Euler implícito===
 
{{matlab|codigo=
 
{{matlab|codigo=
 +
%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)
 
}}
 
}}
 +
 +
[[Archivo:implicitobueno.png|marco|centro|Distribución de la temperatura en el disco]]
 +
 
===Método de Euler modificado===
 
===Método de Euler modificado===
 
{{matlab|codigo=
 
{{matlab|codigo=
 +
%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)
 
}}
 
}}
  
 +
[[Archivo:modificadobueno.png|marco|centro|Distribución de la temperatura en el disco]]
 +
 +
==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 dada por  <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.
 +
{{matlab|codigo=
 +
plot(p,sol(end,:))
 +
}}
 +
 +
[[Archivo:logaritmica.png|marco|centro|Distribución de la temperatura con el tiempo]]
 +
 +
Este es el código para comparar la diferencia con la solución estacionaria
 +
 +
{{matlab|codigo=
 +
%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 off
 +
 +
}}
 +
[[Archivo:1placa.png|marco|centro|Distancia con la solución estacionaria]]
 +
 +
Se observa que la 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.
 +
 +
Para t=0 tenemos:
 +
 +
{{matlab|codigo=
 +
%Datos del problema
 +
clear all
 +
p0=1;
 +
pf=6;
 +
T0=0;
 +
Tf=0;
 +
%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
 +
}}
 +
 +
Para t=1:
 +
{{matlab|codigo=
 +
%Datos del problema
 +
clear all
 +
p0=1;
 +
pf=6;
 +
T0=1;
 +
Tf=1;
 +
%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
 +
}}
 +
Para t=2:
 +
{{matlab|codigo=
 +
%Datos del problema
 +
clear all
 +
p0=1;
 +
pf=6;
 +
T0=2;
 +
Tf=2;
 +
%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
 +
}}
 +
 +
Para t=10:
 +
{{matlab|codigo=
 +
%Datos del problema
 +
clear all
 +
p0=1;
 +
pf=6;
 +
T0=10;
 +
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
 +
}}
 +
 +
{|
 +
|-
 +
| [[Archivo:placa00.png|thumb|500px|left|t=0]] || [[Archivo:placa01.png|thumb|500px|left|t=1|500px]]
 +
|}
 +
 +
{|
 +
|-
 +
| [[Archivo:placa02.png|thumb|500px|left|t=2]] || [[Archivo:placa10.png|thumb|500px|left|t=10|500px]]
 +
|}
 +
 +
En el caso de dividir la discretización en p por 10 observamos:
 +
 +
{{matlab|codigo=
 +
%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
 +
}}
 +
 +
[[Archivo:3placa.png|marco|centro|Con la discretización en p una décima parte]]
 +
==Cambio en la condiciones de frontera==
 +
 +
A continuación cambiaremos a unas condiciones tipo Neumann donde en la frontera exterior colocaremos un material aislante. El aislante hace que no haya pérdida de calor en
 +
ese extremo, es decir, el flujo de temperatura en la dirección radial es nulo.
 +
 +
{{matlab|codigo=
 +
L = 5; dr=0.1; N=L/dr; % discretización espacial
 +
r=1+dr:dr:1+L-dr;
 +
r(1)=0; % lo impongo para la primera iteración
 +
i=2:r;
 +
for n=1:length(r)-1 % bucle
 +
r(i+1)= (1/(1+(1/(n*dr))))*(2*r(i)-r(i-1)+(1/(n*dr))*r(i-1));
 +
end
 +
r(length(r))= r(length(r)-1); % condición de contorno tipo newman
 +
U=r; % solución estacionaria
 +
plot(r)
 +
}}
  
 +
[[Archivo:4placa.png|marco|centro|Solución con las nuevas condiciones de frontera]]
 
[[Categoría:Ecuaciones Diferenciales]]
 
[[Categoría:Ecuaciones Diferenciales]]
 
[[Categoría:ED13/14]]
 
[[Categoría:ED13/14]]
 
[[Categoría:Trabajos 2013-14]]
 
[[Categoría:Trabajos 2013-14]]

Revisión actual del 23:53 19 may 2014

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.

1 Introducción

Placa en anillo con p ∈ [1, 6]

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)
Distribución de la temperatura en el disco

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)


Distribución de la temperatura para p=3

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)


Distribución de la temperatura en el disco

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)


Distribución de la temperatura en el disco

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)


Distribución de la temperatura en el disco

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 dada por [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,:))


Distribución de la temperatura con el tiempo

Este es el código para comparar la diferencia con la solución estacionaria 

%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 off
Distancia con la solución estacionaria

Se observa que la 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.

Para t=0 tenemos: 

%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=0;
%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


Para t=1: 

%Datos del problema
clear all
p0=1;
pf=6;
T0=1;
Tf=1;
%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

Para t=2: 

%Datos del problema
clear all
p0=1;
pf=6;
T0=2;
Tf=2;
%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


Para t=10: 

%Datos del problema
clear all
p0=1;
pf=6;
T0=10;
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


t=0
t=1
t=2
t=10

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


Con la discretización en p una décima parte

4 Cambio en la condiciones de frontera

A continuación cambiaremos a unas condiciones tipo Neumann donde en la frontera exterior colocaremos un material aislante. El aislante hace que no haya pérdida de calor en ese extremo, es decir, el flujo de temperatura en la dirección radial es nulo. 

L = 5; dr=0.1; N=L/dr; % discretización espacial
r=1+dr:dr:1+L-dr;
r(1)=0; % lo impongo para la primera iteración
i=2:r;
for n=1:length(r)-1 % bucle
r(i+1)= (1/(1+(1/(n*dr))))*(2*r(i)-r(i-1)+(1/(n*dr))*r(i-1));
end
r(length(r))= r(length(r)-1); % condición de contorno tipo newman
U=r; % solución estacionaria
plot(r)


Solución con las nuevas condiciones de frontera
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Placa_en_forma_de_Anillo&oldid=13310»