Diferencia entre revisiones de «Ecuacion de calor (Grupo 25C)»

Revisión del 03:46 15 may 2015

1 Introducción

Se considera una placa en forma de anillo comprendida entre los radios ρ=1 y ρ=3 y unos objetos colocados en estas fronteras que mantienen la placa a una temperatura constante de 3 y 0 grados respectivamente.

Representación gráfica de la placa

Consideramos también que su temperatura inicial sigue la siguiente función:

centro

Todas estas características llevan al desarrollo de un problema de calor cuyo sistema se va a plantear a continuación.

2 Sistema de ecuaciones en coordenadas polares

En este caso se va a suponer que la temperatura "U" de la placa solamente depende de la coordenada radial y el tiempo, descartándose así la coordenada cíclica theta. Además se sabe que "U" satisface la siguiente ecuación:

centro

Sabiendo esto se puede desarrollar un sistema de ecuaciones en coordenadas polares que representa el problema planteado.

Al no tener dependencia de la variable “θ” el problema que nos queda es de contorno en “ρ”:

[math] \Delta U = U_{ρρ} +\frac{1 }{ρ} U_ρ [/math]

Aplicando la ecuación que nuestro problema tiene que satisfacer:

[math] U_t - \Delta U = 0 [/math]

[math] U_t - U_{ρρ} -\frac{1 }{ρ} U_ρ = 0 [/math]

Añadiendo las condiciones de contorno y las iniciales nos queda el siguiente planteamiento:

[math] \left\{\begin{matrix} U_t - U_{ρρ} -\frac{1 }{ρ} U_ρ = 0 ρ ∈ (1, 3), t\gt0 \\ U(1,t)=3 U(3,t)=0 t\gt0\\U(ρ,0)=exp(-8*(ρ-2)^2)+3-ρ ρ ∈ (1,3)\end{matrix}\right.[/math]

3 Resolución del problema

Se resuelve por el método de diferencias finitas discretizando la variable espacial "ρ" y la variable temporal "t". En este caso tenemo dos matrices K1 y K2 que acompañan el vector "U" ya que el planteamiento lo exige por las apariciones de las derivadas primera y segunda de U con respecto de "ρ".

3.1 Método del trapecio

En este caso se pide realizar el metodo con los siguientes pasos:

clear all
clc
%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;
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+(k/2)*(K1+K2))\((eye(size(K1))-(k/2)*(K1+K2))*sol(:,i)+k*F');
end
%Grafico
UA=Ua*ones(1,length(t));
UB=Ub*ones(1,length(t));
sol=[UA;sol;UB];
[Mr,Mt]=meshgrid(r,t);
s=surf(Mr,Mt,sol');
set(s,'edgecolor','flat');
colormap summer

Gráfico del problema mediante el método del trapecio

3.2 Temperatura para ρ=2

Para comprobar como cambia la temperatura en función del tiempo para ρ=2 hay que implementar el siguiente codigo en Matlab, a continuación del anterior:

d=(2-a)/h+1;
plot(t,sol(d,:));

Como resultado nos sale el siguiente gráfico donde se observa que para un radio 2, la temperatura baja desde los dos grados hasta los 1.1 para luego mantenerse constante:

Gráfico para ρ=2

3.3 Método de Euler explícito

Para obtener tanto este método como el de Euler implícito y el de Heun, lo unico que hay que modificar es el apartado "%Solucion" del código del método del trapecio.

%Solucion
for i=1:M
    sol(:,i+1)=sol(:,i)+k*(-(K1+K2)*sol(:,i)+F');
end

En este caso se han tenido problemas con la obtención total de los datos y del gráfico para un tiempo de 10 segundos ya que el método es muy impreciso y da valores anormalmente altos. Para un t=1 se ha podido obtener el siguiente gráfico:

Gráfico con el método de Euler explítcito

3.4 Método de Euler implícito

%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end

La gráfica obtenida es muy parecida a la del método del trapecio:

Gráfico meiante el método de Euler implícito

3.5 Método de Heun

El código empleado ha sido el siguiente:

%Solucion
for i=1:M
    Q1=(-(K1+K2))*sol(:,1)+F';
    Q2=(-(K1+K2))*(sol(:,1)+k*Q1)+F';
    sol(:,i+1)=sol(:,1)+(k/2)*(Q1+Q2);
end

Gráfico mediante el método de Heun

4 Planteamiento del problema para tiempos grandes

En este caso es el mismo planteamiento que el anterior pero sin la variable dependiente del tiempo. El problema se reduce a un problema de contorno que se resuelve mediante el siguiente código:

clear all
clc
%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 gráfica resultante nos muestra a qué debería parecerse la función de temperatura para tiempos grandes:

centro

Para t=0

Para t=1

5 Planteamiento del problema con una pieza aislante en ρ=1

6 Consideramos un disco de ρ<3