//  Yves DEBARD - 10 janvier 2010

//  Institut Universitaire de Technologie du MANS
//  Département Génie Mécanique et Productique
//  Avenue Olivier Messiaen
//  72085 LE MANS CEDEX 9

//  Scilab 5.2.0

// ======================================================

// thermique plan : conductivité variable + convection + rayonnement
// (programme non_lin3a)

Stefan=5.67e-8;
TCelsius=273.15;

function [lambda]=Lambda(T)  // T en degrés Celsius
  lambda=5+15*T/500;
endfunction;

L=0.5; // épaisseur en m
// rayonnement sur la face 1
emissivite=0.8;
Tinfini=500+TCelsius;
// convection sur la face 3
h=100;
Tfluide=20+TCelsius;

// vecteur T initial en K
T=[Tinfini;(Tinfini+Tfluide)/2;Tfluide]

// substitution + Newton-Raphson
max_iterations=10;
eps=1e-2;
convergence=%f;
printf('\nItération   dT1     dT2     dT3      ');
printf('T1      T2      T3  (degrés Celcius)\n');
TC=T-TCelsius // vecteur T initial (degrés Celsius)
printf('%5d %35.2f %7.2f %7.2f\n',0,TC(1),TC(2),TC(3));
for i=1:max_iterations
  // matrice de conductivité
  c1=Lambda((TC(1)+TC(2))/2)*2/L;
  c2=Lambda((TC(2)+TC(3))/2)*2/L;
  K=[c1,-c1,0;-c1,c1+c2,-c2;0,-c2,c2];
  // matrice tangente
  Kt=K;
  Kt(1,1)=Kt(1,1)+4*emissivite*Stefan*T(1)^3;
  Kt(3,3)=Kt(3,3)+h;
  // vecteur flux et vecteur résidu
  Flux=[emissivite*Stefan*(Tinfini^4-T(1)^4);0;h*(Tfluide-T(3))]
  Residu=Flux-K*T
  dT=Kt\Residu   // résolution du système linéaire par la méthode de Gauss
  T=T+dT         // mise à jour
  TC=T-TCelsius; // températures en degrés Celsius
  printf('%5d %10.3f %7.3f %7.3f',i,dT(1),dT(2),dT(3));
  printf('%9.2f %7.2f %7.2f\n',TC(1),TC(2),TC(3));
  if norm(dT,'inf')<eps 
    convergence=%t;
    printf('\nConvergence en %d itérations\n\n',i);break;
  end
  if i==max_iterations printf('\nNon convergence en %2d itérations\n\n',i);end
end;

if convergence
  printf('Résidu : ');
  printf('R1 =%6.2f R2 =%6.2f R3 =%6.2f (W/m2)\n',Residu(1),Residu(2),Residu(3));
  printf('Flux : F1 = %5.2f , F3 = %5.2f (W/m2) \n',Flux(1),Flux(3));
  printf('T(2)-0.5*(T(1)+T(3)) = %5.2f degrés Celsius\n',T(2)-0.5*(T(1)+T(3)));
end