%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

%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 = 0.9;     % Emissivité sur la surface coté milieu 1
eps2 = 0.9;     % Emissivité sur la surface coté milieu 2
alpha1 = 0.5 ;  % Absorptivité du matériau coté milieu 1
alpha2 = 0.5 ;  % Absorptivité du matériau coté milieu 2
tho1 = 0.5 ;    % Transmitivité du matériau coté milieu 1
tho2 = 0.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 = matrix(4, 1, 0);

% Initialisation du tableau de sauvegarde
Tsauv = matrix( 4, 1, 1) * 273;
Tnv = matrix( 4, npas, 1) * 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 = 1e4;
while convergence > 1e-3  
%while  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 = matrix( 4, npas, 1) * 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) - Tsauv(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;
        Tsauv(k) = Tnv(k,npas);
    end    
end