clear all
close all
clc

N = 20; %numero de elementos na vertical
M = 30; %numero de elementos na direção radial
niteracoes = 40; %numero de iteracoes
T_agua = zeros (N,niteracoes); %temperatura da agua
T_amb = 273+25; %temperatura ambiente
tempo = zeros(niteracoes,1);
vazao_agua = 0.000025; %[m³/s]
r1 = 2; %raio interno 2cm
rM = 5; %raio externo 5cm
altura = .5; %altura total do cilindro [m]
dz = altura/(N-1); %altura de cada elemento
dt = 1000; %incremento temporal [s]
dr = ((rM - r1)/(M-1)); %tamanho dr de cada elemento (anel)
rho_agua = 1000; %densidade da água
rho = 2700; %aluminio
c_agua = 4200;%calor especifico da agua
c = 920; %calor especifico do material que compoe o cilindro - Aluminio
hc_agua = 1500; %entre 500 e 10000
hc_ar = 50; %entre 10 e 100
k = 237; %para aluminio
sigma = 5.670374419e-8;
epsilon = 0.2;
gerav =0; %10000; %[W/m^3]
IGS_limite = 30;
Tagua_entrada = 273+25;
T = zeros(M,N,niteracoes); %pré alocando espaço
dV = zeros(M,N); %pré alocando espaço
gera = zeros(M,N); %pré alocando espaço
A = zeros(1,M); %pré alocando espaço
qconvar = zeros(N,niteracoes);
qconvagua = zeros(N,niteracoes);
qrad = zeros(N,niteracoes);


%definindo as temperaturas iniciais da agua, superficie interna e miolo
for iz = 1:N
    T_agua(iz,1) = 273+25; %temperatura do elemento iz, no instante inicial
    for P = 1:M
        T(P,iz,1) = 273+25; %temperatura do cilindro
    end

end

%definindo a variação do raio entre os pontos
R(1) = r1;
for P = 2:M
R(P) = R(P-1) + dr;
end

%definindo a variação da area entre os elementos (no caso do anel)
Ainterna = pi*r1*r1;
A(1) = pi*(r1+dr/2)*(r1+dr/2) - Ainterna;
for P = 2:M-1
A(P) = pi*(R(P)+dr/2)*(R(P)+dr/2) - pi*(R(P-1) + dr/2)*(R(P-1) + dr/2);
end
A(M) = pi*(R(M)*R(M)) - pi*(R(M-1) + dr/2)*(R(M-1) + dr/2);

E = 1;
for Q = 1:M
    dV(Q,E) = A(Q)*dz/2;
    gera(Q,E) = dV(Q,E)*gerav;
end

E = N;
for Q = 1:M
    dV(Q,E) = A(Q)*dz/2;
    gera(Q,E) = dV(Q,E)*gerav;
end

for Q = 1:M
for E = 2:N-1
dV(Q,E) = A(Q)*dz;
gera(Q,E) = dV(Q,E)*gerav;
end
end

% Processo iterativo no tempo
for it = 2:niteracoes
    for j=1:N
        T_agua(j,it) = T_agua(j,it-1);
        qconvagua(j,it) = qconvagua(j,it-1);
        qconvar(j,it) = qconvar(j,it-1);
        qrad(j,it) = qrad(j,it-1);
        for p=1:M
        T(p,j,it) = T(p,j,it-1);
        end
  end

    for IGS=1:IGS_limite

        j = 1;
        T_agua(j,it) = (hc_agua*2*pi*r1*(dz/2)*T(1,j,it)+ rho_agua*c_agua*vazao_agua*Tagua_entrada + rho_agua*c_agua*pi*r1*r1*(dz/2)*T_agua(j,it-1)/(dt)) ...
                        /(hc_agua*2*pi*r1*(dz/2) + rho_agua*c_agua*vazao_agua + rho_agua*c_agua*pi*r1*r1*(dz/2)/(dt));



        for j=2:N-1

        T_agua(j,it) = (hc_agua*2*pi*r1*dz*T(1,j,it)+ rho_agua*c_agua*vazao_agua*T_agua(j-1,it) + rho_agua*c_agua*pi*r1*r1*dz*T_agua(j,it-1)/(dt)) ...
                        /(hc_agua*2*pi*r1*dz + rho_agua*c_agua*vazao_agua + rho_agua*c_agua*pi*r1*r1*dz/(dt));

        end

        j = N;
        T_agua(j,it) = (hc_agua*2*pi*r1*(dz/2)*T(1,j,it)+ rho_agua*c_agua*vazao_agua*T_agua(j-1,it) + rho_agua*c_agua*pi*r1*r1*(dz/2)*T_agua(j,it-1)/(dt)) ...
                        /(hc_agua*2*pi*r1*(dz/2) + rho_agua*c_agua*vazao_agua + rho_agua*c_agua*pi*r1*r1*(dz/2)/(dt));

        %para a regiao do "miolo" da capa cilindrica (conducao apenas)

        for P=2:M-1
            for j=2:N-1
                T(P,j,it) = ((2*pi*dz*k/log(R(P)/R(P-1)))*T(P-1,j,it) + (2*pi*dz*k/log(R(P+1)/R(P)))*T(P+1,j,it)+ (k*A(P)/dz)*T(P,j-1,it)+ (k*A(P)/dz)*T(P,j+1,it) + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*dz*k/log(R(P)/R(P-1))) +(2*pi*dz*k/log(R(P+1)/R(P)))+ k*A(P)/dz + k*A(P)/dz + rho*c*dV(P,j)/dt);
            end
    end

        %para a regiao lateral externa (conducao, radiacao e conveccao com o ar)
        P=M;
            for j=2:N-1
                C_r(j) = sigma*epsilon*2*pi*R(P)*dz*(T(P,j,it)^2 + T_amb^2)*(T(P,j,it)+T_amb);
                T(P,j,it) = ((2*pi*dz*k/log(R(P)/R(P-1)))*T(P-1,j,it) + (k*A(P)/dz)*T(P,j-1,it)+(k*A(P)/dz)*T(P,j+1,it) + hc_ar*2*pi*R(P)*dz*T_amb + C_r(j)*T_amb + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*dz*k/log(R(P)/R(P-1))) + k*A(P)/dz + k*A(P)/dz + hc_ar*2*pi*R(P)*dz + C_r(j) + rho*c*dV(P,j)/dt);
                qconvar(j,it) = hc_ar*2*pi*R(P)*dz*(T_amb-T(P,j,it));
                qrad(j,it) = C_r(j)*(T_amb-T(P,j,it));
            end


        %para a regiao lateral interna (conducao e conveccao com a agua)
        P=1;
            for j=2:N-1
                T(P,j,it) = ((2*pi*dz*k/log(R(P+1)/R(P)))*T(P+1,j,it) + (k*A(P)/dz)*T(P,j-1,it)+(k*A(P)/dz)*T(P,j+1,it) + hc_agua*2*pi*R(P)*dz*T_agua(j,it) + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*dz*k/log(R(P+1)/R(P))) + k*A(P)/dz + k*A(P)/dz + hc_agua*2*pi*R(P)*dz + rho*c*dV(P,j)/dt);
                       qconvagua(j,it) = hc_agua*2*pi*r1*(dz)*(T(1,j,it)-T_agua(j,it));
            end


        %para a regiao superior (situacao adiabatica superficie de cima)
    for P=2:M-1
            for j=1
                T(P,j,it) = ((2*pi*(dz/2)*k/log(R(P)/R(P-1)))*T(P-1,j,it) + (2*pi*(dz/2)*k/log(R(P+1)/R(P)))*T(P+1,j,it) + (k*A(P)/dz)*T(P,j+1,it) + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*(dz/2)*k/log(R(P)/R(P-1))) +(2*pi*(dz/2)*k/log(R(P+1)/R(P))) + k*A(P)/dz + rho*c*dV(P,j)/dt);
            end
    end

         %para a regiao inferior (situacao adiabatica superficie de baixo)
    for P=2:M-1
            for j=N
                T(P,j,it) = ((2*pi*(dz/2)*k/log(R(P)/R(P-1)))*T(P-1,j,it) + (2*pi*(dz/2)*k/log(R(P+1)/R(P)))*T(P+1,j,it) + (k*A(P)/dz)*T(P,j-1,it) + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*(dz/2)*k/log(R(P)/R(P-1))) +(2*pi*(dz/2)*k/log(R(P+1)/R(P))) + k*A(P)/dz + rho*c*dV(P,j)/dt);
            end
    end

        %para a regiao superior + canto esquerdo (situacao adiabatica e contato com agua)
        P=1;
                j=1;
                T(P,j,it) = ((2*pi*(dz/2)*k/log(R(P+1)/R(P)))*T(P+1,j,it) + (k*A(P)/dz)*T(P,j+1,it) + hc_agua*2*pi*R(P)*(dz/2)*T_agua(j,it) + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*(dz/2)*k/log(R(P+1)/R(P))) + k*A(P)/dz + hc_agua*2*pi*R(P)*(dz/2) + rho*c*dV(P,j)/dt);
                 qconvagua(j,it) = hc_agua*2*pi*r1*(dz/2)*(T(1,j,it)-T_agua(j,it));



        %para a regiao superior + canto direito (situacao adiabatica e contato com ar)
        P=M;
                j=1;
                C_r(j) = sigma*epsilon*2*pi*R(P)*(dz/2)*(T(P,j,it)^2 + T_amb^2)*(T(P,j,it)+T_amb);
                T(P,j,it) = ((2*pi*(dz/2)*k/log(R(P)/R(P-1)))*T(P-1,j,it) + (k*A(P)/dz)*T(P,j+1,it) + hc_ar*2*pi*R(P)*(dz/2)*T_amb + C_r(j)*T_amb + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*(dz/2)*k/log(R(P)/R(P-1))) + k*A(P)/dz + hc_ar*2*pi*R(P)*(dz/2) + C_r(j) + rho*c*dV(P,j)/dt);
                qconvar(j,it) = hc_ar*2*pi*R(P)*(dz/2)*(T_amb-T(P,j,it));
                qrad(j,it) = C_r(j)*(T_amb-T(P,j,it));



        %para a regiao inferior + canto esquerdo (situacao adiabatica e contato com agua)
        P=1;
                j=N;
                T(P,j,it) = ((2*pi*(dz/2)*k/log(R(P+1)/R(P)))*T(P+1,j,it) + (k*A(P)/dz)*T(P,j-1,it) + hc_agua*2*pi*R(P)*(dz/2)*T_agua(j,it) + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*(dz/2)*k/log(R(P+1)/R(P))) + k*A(P)/dz + hc_agua*2*pi*R(P)*(dz/2) + rho*c*dV(P,j)/dt);
                  qconvagua(j,it) = hc_agua*2*pi*r1*(dz/2)*(T(1,j,it)-T_agua(j,it));


        %para a regiao inferior + canto direito (situacao adiabatica e contato com ar)
        P=M;
                j=N;
                C_r(j) = sigma*epsilon*2*pi*R(P)*(dz/2)*(T(P,j,it)^2 + T_amb^2)*(T(P,j,it)+T_amb);
                T(P,j,it) = ((2*pi*(dz/2)*k/log(R(P)/R(P-1)))*T(P-1,j,it) + (k*A(P)/dz)*T(P,j-1,it) + hc_ar*2*pi*R(P)*(dz/2)*T_amb + C_r(j)*T_amb + gera(P,j) + (rho*c*dV(P,j)*T(P,j,it-1)/dt)) ...
                            /((2*pi*(dz/2)*k/log(R(P)/R(P-1))) + k*A(P)/dz + hc_ar*2*pi*R(P)*(dz/2) + C_r(j) + rho*c*dV(P,j)/dt);
                qconvar(j,it) = hc_ar*2*pi*R(P)*(dz/2)*(T_amb-T(P,j,it));
                qrad(j,it) = C_r(j)*(T_amb-T(P,j,it));


    end

    tempo(it)=tempo(it-1)+dt;
end