%programme de bilan enthalpique sur un composant discrétisé en 4 noeuds
%interne (dont deux de surface).
%Les conditions limites sont définies dans les milieux 1 et 2;
%
%Le schéma equivalent est le suivant
clear all;
clc;

tic;

%Apports radiatifs  (IR1)        |                         | IR2
%Echanges radiatifs (Tr1 et hr1) | T1-----T2-----T3-----T4 | Tr2 et hr2
%Ecanges convectifs (Te1 et hc1) |                         | Te2 et hc2

%Définition des conditions limites
%-------------milieu 1
IR1 = 0;       % apports radiatifs (W/m²)
Tr1 = 30+273;      % Température equivalente pour les échanges radiatifs (°C)
Te1 = 30+273;      % Température équivalente pour les échange convectifs (°C)
hc1 = 20;      % Coefficient d'échange convectif coté mimieu 1 (W/m²/K)

%-------------milieu 2
IR2 = 0;       % apports radiatifs (W/m²)
Tr2 = 10+273;      % Température equivalente pour les échanges radiatifs (°C)
Te2 = 10+273;      % Température équivalente pour les échange convectifs (°C)
hc2 = 20;      % Coefficient d'échange convectif coté mimieu 2 (W/m²/K)

%Définition des caractéristiques radiatives
eps1 = .9;     % Emissivité sur la surface coté milieu 1
eps2 = .9;     % Emissivité sur la surface coté milieu 2
alpha1 = .5 ;  % Absorptivité du matériau coté milieu 1
alpha2 = .5 ;  % Absorptivité du matériau coté milieu 2
tho1 = .5 ;    % Transmitivité du matériau coté milieu 1
tho2 = .5 ;    % Transmitivité du matériau coté milieu 2

%Définition des caractiristiques du matériau
lambda = 0.201;  % Conductivité thermique (W/m/K)
rho = 400;      % densité (kg/m3)
Cp = 1400;       % Capacité calorifique (J/kg/K)
ep = 5E-2;       % Epaisseur du matériau (m)


%Définition des conditions initiales
Tinit = 0;       % Température de condition initale (°C)

%definition des cond et capacité en chaque noeud
cond = lambda / ep * 3;
capa = rho * Cp * ep / 3;
capa2 = capa / 2;
condneg = -1 * cond;  %matrice négative des conductances
boltz = 5.67e-8;

% Temps final physique
tempsfinal = 1 * 3600;  

% Données numériques
npas = 1;
NNpas = 1E+3;
dt = tempsfinal / npas;
conv = zeros( 4, 1 );

% Initialisation du tableau de sauvegarde
T = ones( 4, 1) * inf;
Tnv = ones( 4, npas ) * Tinit + 273;

% Choix de la progression 
    % Geométrique = 0;
    % Arithmétique != 0;
rgeo = 2;  % Raison de la progression géométrique
rarithm = 1;  % Raison de la progression arithmétique
progression = 1;


%programme de convergence des températures
convergence = inf;
while ( ( convergence > 1E-3 ) && ( npas < NNpas ) )
    
    % Sélection de la progression
    if ( progression )
        % Incrémentation du nombre de pas de temps
        npas = npas + rarithm;
    else
        npas = npas * rgeo;
    end
    
    % Calcul de la valeur du pas de temps
    dt = tempsfinal / ( npas - 1 );
    
    % Initialisation
    Tnv = ones( 4, npas ) * Tinit + 273;
    
    for n = 2 : 1 : npas

        %calcul des coefficients d'échange radiatifs
        hr1 = eps1 * boltz * ( Te1^2 + Tnv(1,n-1)^2 ) * ( Te1 + Tnv(1,n-1) );
        hr2 = eps2 * boltz * ( Te2^2 + Tnv(4,n-1)^2 ) * ( Te2 + Tnv(4,n-1) );

        %definition de la matrice des capacités
        MAT(1,1) = capa2 / dt + hr1 + hc1 + cond;
        MAT(1,2) = condneg;
        MAT(1,3) = 0;
        MAT(1,4) = 0;
        MAT(2,1) = condneg;
        MAT(2,2) = capa / dt + 2 * cond;
        MAT(2,3) = condneg;
        MAT(2,4) = 0;
        MAT(3,1) = 0;
        MAT(3,2) = condneg;
        MAT(3,3) = capa / dt + 2 * cond;
        MAT(3,4) = condneg;
        MAT(4,1) = 0;
        MAT(4,2) = 0;
        MAT(4,3) = condneg;
        MAT(4,4) = capa2/ dt + cond + hr2 + hc2;

        % Définition du veteur résultat
        SOL(1) = capa2 / dt * Tnv(1,n-1) + hc1 * Te1 + hr1 * Tr1 + alpha1 * IR1;
        SOL(2) = capa / dt * Tnv(2,n-1);
        SOL(3) = capa / dt * Tnv(3,n-1);
        SOL(4) = capa2 / dt* Tnv(4,n-1) + hc2 * Te2 + hr2 * Tr2 + alpha2 * IR2;

        % Résolution du système
        IMAT = inv(MAT);
        SOLt = [SOL(1),SOL(2),SOL(3),SOL(4)];
        for k = 1 : 1 : 4;
        Tnv(k,n) = IMAT(k,1) * SOLt(1) + IMAT(k,2) * SOLt(2) + IMAT(k,3) * SOLt(3) + IMAT(k,4) * SOLt(4);
        end    
    end
    
    % Calcul de la norme infinie
    for k = 1 : 1 : 4;
        conv(k) = Tnv(k,npas) - T(k);
    end
    
    % Calcul du carré de la norme euclidienne d'ordre 2
    convergence = 0.;
    for n = 1 : 1 : 4
        convergence = convergence + conv(n)^2;
    end
    
    % Sauvegarde du champ de température
    for k = 1 : 1 : 4;
        T(k) = Tnv(k,npas);
    end    
end

duree = toc;
fprintf( '\n\n Temps de calcul : %.2E s', duree );

% Abscisse de traçage des graphiques
temps = linspace( 0, tempsfinal, npas );

% Ordonnées initiales
VTe1 = ones(1,length(Tnv(4,:))).*Te1-273;
VTe2 = ones(1,length(Tnv(4,:))).*Te2-273;

Tnv = Tnv - 273;
figure(1);
clf;
hold on;
plot( temps, Tnv(1,:),'g', temps,Tnv(2,:),'y',temps,Tnv(3,:),'b',temps,Tnv(4,:),'r')
plot( temps, VTe1, 'k','linewidth', 2 );
plot( temps, VTe2, 'k', 'linewidth', 2 );
grid on;
box on;
xlabel('Temps (h)');
ylabel('Température (°C)');
legend('Tsurf1','Tint12','Tint21','Tsurf2','Tmilieu1','Tmilieu2', 'location', 'best' );
axis( [ 0 tempsfinal min(Tinit,5) max(Te1,Te2)-273+5 ] );

fprintf( '\n\n\t\tFIN\n\n' );