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 2 t0=0; 3 tM=10; %instante final 4 q=input(' Introduce difusividad'); %1 % difusividad= conductividad termica/densidad y calor especifico 5 6 g=@(r,t) 0*r; 7 ca=@(t) 3+0*t; 8 cb=@(t) 0*t; 9 u0=@(r) exp(-8*(r-2).^2)+3-r; 10 11 % discretizacion espacial 12 h=0.1;%Anchura de paso en el espacio 13 N=round((3-1)/h); %subintervalos en x 14 r=linspace(1,3,N+1); 15 r=r'; % El vector r tiene N+ 1 elementos 16 17 % Ejerciocio tipo Condiciones Dirichlet 18 19 rr=r(2:end-1); % tomamos los elementos del segundo al penultimo (N-1 ecuaciones)(N componentes) 20 21 %Matriz K de N-1 x N-1) 22 K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); 23 24 L=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1)); 25 A=diag(1./rr); 26 L=A*L; 27 28 %Valor inicial 29 U0=u0(rr); 30 31 % Vector de terminos independientes 32 % depende del tiempo ---> queda para despues 33 34 k=h/4; 35 M=round((tM-t0)/k); % numero de subintervalos en t 36 t=linspace(t0,tM,M+1); % vector con los tiempos 37 38 % sol matriz que tendra la solucion para cada punto e intante 39 sol=zeros(N-1,M+1); 40 sol(:,1)=U0; 41 42 % Aplicamos el metodo de Trapecio 43 for i=1:M 44 G1=g(rr,t(i)); 45 G1(1)=G1(1)+q*ca(t(i))/h^2-(ca(t(i)))/(2*h*rr(1)); 46 G1(end)=G1(end)+q*cb(t(i))/h^2+(cb(t(i)))/(2*h*rr(end)); 47 G2=g(rr,t(i+1)); 48 G2(1)=G2(1)+q*ca(t(i+1))/h^2-(ca(t(i)))/(2*h*rr(1)); 49 G2(end)=G2(end)+q*cb(t(i+1))/h^2+(cb(t(i)))/(2*h*rr(end)); 50 sol(:,i+1)=(eye(size(K))+(k/2)*(K+L))\(sol(:,i)+k/2*(-(K+L)*sol(:,i)+G1+G2)); 51 end 52 53 % tenemos la solucion completa SALVO las condiciones DIRICHLET 54 Ua=ca(t); 55 Ub=cb(t); 56 sol=[Ua;sol;Ub]; % las añadimos 57 58 % gráfico 59 60 [Mt,Mr]=meshgrid(t,r); 61 62 surf(Mr,Mt,sol) 63 64 xlabel('x'); 65 ylabel('t');
