% ALETTA
% Calcolo della distribuzione di temperatura lungo una 
% aletta circolare con il metodo delle differenze finite
% e valutazione della efficienza e della efficacia della aletta


d=input('Diametro aletta? (m)_');
k = input('Conducibilità termica? (W/m.K)_');
L = input('Lunghezza aletta? (m)_'); 
h = input('Coeff. di scambio termico superficiale? (W/m2.K)_');
n = input('Numero di suddivisioni? _');

dx = L/n; % incremento deltax
Kk = k*pi*(d^2)/(4*dx); % trasmittanza conduttiva
Kc = h*pi*d*dx; % trasmittanza convettiva
Kce = h*pi*(d^2)/4; % trasmittanza convettiva di estremità
OK = 0;
chk = 0.01; % approssimazione

Tb = input('Temperatura della base? (°C)_');
Tinf = input('Temperatura del fluido? (°C)_');

T0=Tb-Tinf; %Temperatura della base relativizzata rispetto a quella del fluido
T(1,1) = T0; % Temperatura del primo nodo FISSA

T(2:n+1,1) = T0; % Assegnazione delle temperature iniziali

while OK == 0
   i = 2:n;
   % Calcolo delle temperature nodali T(i,2)in base a quelle assegnate
   % o calcolate nello step precedente e contenute in T(i,1)
   T(i,2) = (T(i-1,1)+T(i+1,1))/(2+(Kc/Kk)); % Calcolo delle temperature nodali
   T(1,2) = T0; % Temperatura di base FISSA
   T(n+1,2) = T(n,1)/(1+Kc/(2*Kk)+Kce/Kk); % Temperatura del nodo sommitale
   
   DT = abs(T(:,1)-T(:,2)); % Calcolo della differenza fra le temperature di due step consecutivi
  
   % Trasferimento dei valori di T appena calcolati nel vettore
   % delle temperature di riferimento (step precedente)
   T(:,1)=T(:,2); 
   
   %Verifica della differenza fra le temperature
   if DT < chk %Se è minore di quella richiesta ESCI
      OK = 1;
   end
end

Tf=T(:,1)+Tinf % Distribuzione finale di temperatura lungo l'aletta
Clear T

% Routine di plottaggio
plot((0:n)*L/n,Tf(1:n+1))
axis tight
grid on
xlabel('Lunghezza (m)')
ylabel('Temperatura (°C)')
%Calcolo della potenza termica dissipata attraverso l'aletta
Qx(1) = (Tf(1)-Tinf)*Kc/2;
Qx(2:n) = (Tf(n)-Tinf)*Kc;
Qx(n+1) = (Tf(n+1)-Tinf)*(Kc/2+Kce);
Qa = sum(Qx(1:n+1));

% Potenza termica che l'aletta dissiperebbe se si trovasse a Tb costante
Qt = (Tb-Tinf)*(pi*d*L+pi*(d^2)/4)*h;

efficienza=Qa/Qt


Qb = (Tb-Tinf)*h*pi*(d^2)/4; % Potenza termica dissipata senza aletta
Qba = (Tb-Tf(2))*Kk + (Tb - Tinf)*Kc/2; % Potenza termica dissipata con aletta 
efficacia = Qba/Qb

clear all