* fichier :  echang.dgibi

************************************************************************
*                                                                      *
* Echange par convection bilatéral                                     *
*                                                                      *
* Calcul 3D (ou 2D-plan) linéaire permanent                                         *
*                                                                      *
* on considère 3 plaques parallèles                                    *
*                                                                      *
*    -------- T3                                                       *
*                      face à face                                     *
*                                                                      *
*    -------- T2                                                       *
*                      face à face                                     *
*                                                                      *
*    -------- T1,e1                                                    *
*                                                                      *
* La température des plaques supérieure et inférieure est imposée et   *
* uniforme (respectivement 1000K et 2000K) et on détermine la tempéra- *
* ture d'équilibre de la plaque intermédiaire.                         *
*                                                                      *
************************************************************************

OPTI ECHO 0;

xdim3 = VRAI;


*-------------------------------
* Données physiques du problème
*-------------------------------

T1 = 1000. ; T3 = 2000. ;
h12 = 200. ; h23 = 200. ;
lamb = 2.e3 ;

*-------------------------------
* Solution analytique
*-------------------------------

T2_ana =  ((h12*T1) + (h23*T3))/(h12 + h23) ;

*-------------------------------
*   Maillage 
*------------------------------- 

* 3D

si xdim3 ;
option  dime 3 elem cub8  ;
p1=0 0 0 ; p2=1 0 0 ; p3=1 1 0 ; p4=0 1 0;
n= 3 ;
l1=p1 d n p2;
l2=p2 d n p3;
l3=p3 d n p4;
l4=p4 d n p1;
s1= (l1 et l2 et l3 et l4) surf plan ;

n=1 ;
v = 0 0 0.1 ;  v1 = (-1)*v ;
w = 0 0 0.001 ;  w1 = (-1)*w ;

 s2inf = s1 ;
 s2sup = s2inf plus w ;
 plaq2 = s2inf volu n s2sup ;

 s3inf = s2sup plus v ;
 s3sup = s3inf plus w ;
 plaq3 = s3inf volu n s3sup ;

 s1sup = s2inf plus v1;
 s1inf = s1sup plus w1;
 plaq1 = s1INF volu n s1sup ;

 rac12 = liai  0.1  s1sup s2inf ;
 rac23 = liai  0.1  s2sup s3inf ;

 tout = plaq1 et plaq2 et plaq3 ;


* 2D

sinon ;
  option  dime 2 elem qua4  ;
  p  = 0. 0. ; q  = 1. 0. ; s  = droi 1 p q ;
  s1inf = s ;
  plaq1 = s1inf tran 1 (0. 0.001) ;
  s1sup = inve ( cote 3 plaq1) ;

  s2inf = s1inf plus (0. 0.1) ;
  plaq2 = s2inf tran 1 (0. 0.001) ;
  s2sup = inve ( cote 3 plaq2) ;

  s3inf = s2sup plus (0. 0.1) ;
  plaq3 = s3inf tran 1 (0. 0.001) ;

  rac12 = racc 0.1 s1sup s2inf ;
  rac23 = racc 0.1 s2sup s3inf ;

  tout = plaq1 et plaq2 et plaq3 ;

finsi ;

*------------------------------- 
*  Conduction
*------------------------------- 

  mcd = mode tout thermique ;
  k = mate mcd 'K' lamb ;
  cnd = cond mcd k ;

*------------------------------- 
* Convection bilatérale
*------------------------------- 

  mcv12 = mode rac12 thermique convection ;
  mcv23 = mode rac23 thermique convection ;

  ma12 = matr mcv12 'H' h12 ;
  ma23 = matr mcv23 'H' h23 ;

  cr12 = cond mcv12 ma12    ;
  cr23 = cond mcv23 ma23    ;

*------------------------------- 
*  Conditions aux limites
*------------------------------- 

    c1 = bloq plaq1   'T' ;
    tim1 = depi c1 T1 ;

    c3 = bloq plaq3   'T' ;
    tim3 = depi c3 T3 ;

*------------------------------- 
*  Résolution
*------------------------------- 


       crtot= cr12 et cr23  ;
       cndtot = (crtot et cnd) et (c1 et c3) ;

       tt = resou cndtot (tim3 et tim1) ;

    Temp = exco 'T' tt 'T' ;
    T2 = redu plaq2   Temp ;


*------------------------------- 
*  code fonctionnement
*------------------------------- 
   nno = nbno plaq2;
   rr = resul t2 ;
   rra = MAXI ( rr / nno );
   ere = abs  (( rra - T2_ana ) / T2_ana );
   MESS ' ERREUR RELATIVE '   ere;


   SI (ere  <EG 1.d-8 ) ;
     ERRE 0;
   SINON;
     ERRE 5;
   FINSI;
fin;