Diferencia entre revisiones de «Placa en forma de Anillo»

De MateWiki
Saltar a: navegación, buscar
(→‎Evolución de la solución en el tiempo)
Línea 320: Línea 320:
 
[[Archivo:logaritmica.png|marco|centro|Distribución de la temperatura con el tiempo]]
 
[[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=
 
{{matlab|codigo=
Línea 355: Línea 355:
  
 
}}
 
}}
−
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.
+
[[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:
 
Para t=0 tenemos:
Línea 493: Línea 495:
 
}}
 
}}
  
 +
[[Archivo:2placa.png|marco|centro|Para los diferentes t]]
  
 
En el caso de dividir la discretización en p por 10 observamos:
 
En el caso de dividir la discretización en p por 10 observamos:
Línea 529: Línea 532:
 
}}
 
}}
  
 +
[[Archivo:3placa.png|marco|centro|Con la discretización en p una décima parte]]
 
==Cambio en la condiciones de frontera==
 
==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�erdida de calor en
+
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�on radial es nulo.
+
ese extremo, es decir, el flujo de temperatura en la dirección radial es nulo.
  
 
{{matlab|codigo=
 
{{matlab|codigo=
Línea 546: Línea 550:
 
plot(r)
 
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 del 23:45 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


Para los diferentes t

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=13288»