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 | |
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.
2 Resolución del problema
Para resolver problema se emplea el método de diferencias finitas con la siguiente discretización: Δρ=0.1, Δt=Δρ/4 en t∈[0,10] 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
[1]- ↑ % 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 % Ejerciocio 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');
