% PARETE
% Applicazione per il calcolo della distribuzione di temperatura 
% entro una parete composta, con il metodo delle differenze finite
% I parametri possono essere modificati solo all'interno dello
% script


nx = 11; %numero di righe
ny = 31; %numero di colonne
dx = 0.04; %passo dx (verticale)
dy = 0.04; %passo dy (orizzontale)
Lx=(nx-1)*dx;
Ly=(ny-1)*dy;

OK = 0;
chk = 0.05; % approssimazione

T1 = 20; %Temperatura sulla faccia inferiore
T2 = 0; %Temperatura sulla faccia superiore
h1 = 7.7; %Coeff. di scambio convettivo inferiore
h2 = 25; %Coeff. di scambio convettivo superiore

k1 = 0.7; % mattoni
k2 = 2.5; % cls
k3 = 0.05; % isolante

% Per assegnare le conducibilità ai nodi 
% modificare le righe sottostanti individuando le aree
% omogenee

Kk(1:nx,1:ny) = k1; %tutto omogeneo
Kk(1:nx,11:20) = k2; %striscia intermedia di cls (pilastro)
Kk(10:11,1:ny) = k3; %fascia di isolante superiore
% ----------------------------------------------------

% Plottaggio della mappa delle conducibilità
contourf(Kk(1:nx,1:ny),10)
colorbar
pause

% Calcolo delle trasmittanze relative a ogni nodo *******
% Kz(n,m,i) i da 1 a 4

Kz(1,1:ny,1)=h1*dy; %Resistenza liminare interna
Kz(2:nx,1:ny,1)=(Kk(2:nx,1:ny)+Kk(1:nx-1,1:ny))/2;

Kz(1:nx,1:ny-1,2)=(Kk(1:nx,1:ny-1)+Kk(1:nx,2:ny))/2;
Kz(1:nx,ny,2)=0; % Piano adiabatico
Kz(1,1:ny,2)=Kz(1,1:ny,2)*0.5; %spessore dimezzato
Kz(nx,1:ny,2)=Kz(nx,1:ny,2)*0.5; %spessore dimezzato

Kz(1:nx-1,1:ny,3)=(Kk(1:nx-1,1:ny)+Kk(2:nx,1:ny))/2;
Kz(nx,1:ny,3)=h2*dy; %Resistenza liminare esterna

Kz(1:nx,2:ny,4)=(Kk(1:nx,2:ny)+Kk(1:nx,1:ny-1))/2;
Kz(1:nx,1,4)=0; %Piano adiabatico
Kz(1,1:ny,4)=Kz(1,1:ny,4)*0.5; %spessore dimezzato
Kz(nx,1:ny,4)=Kz(nx,1:ny,4)*0.5; %spessore dimezzato

% assegnazione temperatura iniziale
T(1:nx,1:ny)=(T2+T1)/2;

% routine di calcolo iterativo
while OK == 0
   
   Tz(1,1:ny,1)=T1;
   Tz(2:nx,1:ny,1)=T(1:nx-1,1:ny);
   Tz(1:nx,1:ny-1,2)=T(1:nx,2:ny);
   Tz(1:nx-1,1:ny,3)=T(2:nx,1:ny);
   Tz(nx,1:ny,3)=T2;
   Tz(1:nx,2:ny,4)=T(1:nx,1:ny-1);
   for i=1:nx
      for j=1:ny
         KT(i,j)=sum(Kz(i,j,1:4).*Tz(i,j,1:4));
         Tn(i,j)=KT(i,j)/sum(Kz(i,j,1:4));
      end
   end
   DT = abs(Tn(1:nx,1:ny)-T(1:nx,1:ny));
   T(1:nx,1:ny)=Tn(1:nx,1:ny);
  
   if DT < chk
      OK = 1;
   end
end

% Plottaggio risultati-------------------------
[px,py]=gradient(T(1:nx,1:ny));
px=-px; py=-py;
%pcolor(T(1:nx,1:ny))
contourf(T(1:nx,1:ny),20)
%shading flat
axis tight
colormap('jet')
colorbar
%[c,h]=contourf(T(1:nx,1:ny),10); clabel(c,h,'manual')
hold on

pause
% disegna le linee di flusso
quiver(px,py)
hold off


%Qi(1:ny)=h1*dy*(T1-T(1,1:ny));
%Qt=sum(Qi)