Download couplage_TH1D_Th3D_2.dgibi

Back to the list

* fichier :  couplage_TH1D_Th3D_2.dgibi
************************************************************************
* NOM         : couplage_TH1D_Th3D_2.dgibi
* DESCRIPTION :
*
* Voir rapport DM2S/SFME/LTMF/RT/03-012/B
* Modélisation couplée thermique 3D-thermohydraulique 1D d'un coeur de
* réacteur à caloporteur gaz.
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 29/04/2003, version initiale
* HISTORIQUE : v1, 29/04/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
* interact = VRAI : pour une utilisation interactive
* graph    = VRAI : pour tracer des choses
* irrad    = VRAI : si le combustible est irradié
* pvar     = VRAI : si la puissance dégagée dans le combustible à un
*                   profil sinusoïdale (en fonction de z)
* discr    = LINE ou QUAF : discrétisation linéaire (suffisante)
*                           ou quadratique
interact= FAUX ;
graph   = FAUX ;
irrad   = VRAI ;
pvar    = VRAI ;
discr   = 'LINE' ;
*discr = 'QUAF' ;
************************************************************************
*
* Début des procédures
*
'OPTION' 'ECHO' 0 ;
* NOM         : ERRREL
* DESCRIPTION : Calcul d'une erreur relative
'DEBPROC' ERRREL ;
'ARGUMENT' val*'FLOTTANT' ;
'ARGUMENT' valref*'FLOTTANT' ;
*
'SI' ('<' ('ABS' valref) 1.D-10) ;
   echref = 1.D0 ;
'SINON' ;
   echref = valref ;   
'FINSI' ;
*
errabs = 'ABS' ('/' ('-' val valref) echref);
*
'RESPRO' errabs ;
'FINPROC' ;
* NOM         : TRACCHPO
* DESCRIPTION : Tracé d'un chpoint avec titre optionnel.
'DEBPROC' TRACCHPO ;
'ARGUMENT' pn*'CHPOINT ' ;
'ARGUMENT' mt*'MAILLAGE' ;
'ARGUMENT' mt2/'MAILLAGE' ;
'ARGUMENT' oeil/'POINT   ' ;
'ARGUMENT' tit/'MOT     ' ;
'SI' ('EXISTE' tit) ;
   titpn = 'CHAINE'  'FORMAT' formflot tit titglob ;
'SINON' ;   
   titpn = 'CHAINE' 'FORMAT' formflot titglob ;
'FINSI' ;
rescal = pn ;
'SI' ('EXISTE' oeil) ;
   'SI' ('EXISTE' mt2) ;
      'TRACER' 'CACH' oeil rescal mt mt2 'TITR' titpn nbisov ;
   'SINON' ;   
      'TRACER' 'CACH' oeil rescal mt mt 'TITR' titpn nbisov ;
   'FINSI' ;
'SINON' ;
   'SI' ('EXISTE' mt2) ;
      'TRACER' 'CACH' rescal mt mt2 'TITR' titpn nbisov ;
   'SINON' ;   
      'TRACER' 'CACH' rescal mt mt 'TITR' titpn nbisov ;
   'FINSI' ;
'FINSI' ;
'FINPROC' ;
* NOM         : TRACCHML
* DESCRIPTION : Trace un CHPOINT défini sur les centres des éléments
*               (valeur constante par élément) avec titre optionnel
'DEBPROC' TRACCHML ;
'ARGUMENT' chml*'CHPOINT ' ;
'ARGUMENT' $mt*'MMODEL  ' ;
'ARGUMENT' mt2/'MAILLAGE' ;
'ARGUMENT' oeil/'POINT   ' ;
'ARGUMENT' tit/'MOT     ' ;
'SI' ('EXISTE' tit) ;
   titpn = 'CHAINE'  'FORMAT' formflot tit titglob ;
'SINON' ;   
   titpn = 'CHAINE' 'FORMAT' formflot titglob ;
'FINSI' ;
rescal = 'KCHA' $mt chml 'CHAM' ;
modbid = 'MODELISER' ('DOMA' $mt 'MAILLAGE') 'THERMIQUE' ;
'SI' ('EXISTE' oeil) ;
   'SI' ('EXISTE' mt2) ;
      'TRACER' 'CACH' oeil rescal modbid mt2 'TITR' titpn nbisov ;
   'SINON' ;   
      'TRACER' 'CACH' oeil rescal modbid mt 'TITR' titpn nbisov ;
   'FINSI' ;
'SINON' ;
   'SI' ('EXISTE' mt2) ;
      'TRACER' 'CACH' rescal modbid mt2 'TITR' titpn nbisov ;
   'SINON' ;   
      'TRACER' 'CACH' rescal modbid mt 'TITR' titpn nbisov ;
   'FINSI' ;
'FINSI' ;   
'FINPROC' ;
* NOM         : PRPHE
*                Proprietes thermiques de l Helium 4
* Input : PTOT : pression en Pa ou en bars
*         T ou TS : temperature sous forme de CHPOINT ou FLOTTANT
*
* Output : RHOHE : densite de l helium (USI) de meme type que T ou TS
*          MUHE  : viscosite
*          LHE   : conductivite thermique
*          CPHE  : capacité calorifique a pression constante
'DEBPROC' PRPHE PTOT*'FLOTTANT' T/'CHPOINT' TS/'FLOTTANT' ;
PREF = 101325. ;
'SI' ('<EG' PTOT PREF) ;
     P = PTOT '*' PREF ;
'FINSI' ;
RHE = 2077.22 ;
CPHE = 5193. ;
C1 = 9.489433E-04 ;
C2 = 9.528079E-04 ;
C3 = 3.420680E-02 ;
C4 = 2.739470E-03 ;
C5 = 9.409120E-04 ;
'SI' ('EXIST' T ) ;
     B = 'INVE' (T '*' C5 '+' 1.0) ;
     B = B '*' C4 '-' (( 'INVE' (T '*' C3 '-' 1.0)) '*' C2);
     B = B '+' C1 ;
     RHOHE = ('INVE' (T '*' RHE '+' (B '*' P))) '*' P ;
     MUHE = ( T '**' 0.687) '*' 3.953E-07 ;
     LHE = ( T '**' 0.701 ) '*' 2.774E-03 ;
     'RESPROC' RHOHE MUHE LHE CPHE ;
'FINSI' ; 
'SI' ('EXIST' TS ) ;
     B = 1.0 '/' (TS '*' C5 '+' 1.0) ;
     B = B '*' C4 '-' (C2 '/' (TS '*' C3 '-' 1.0)) ;
     B = B '+' C1 ;
     RHOHE = (P '/'  (TS '*' RHE '+' (B '*' P)));
     MUHE = ( TS '**' 0.687) '*' 3.953E-07 ;
     LHE = ( TS '**' 0.701 ) '*' 2.774E-03 ;
     'RESPROC' RHOHE MUHE LHE CPHE ;
'FINSI' ; 
'FINPROC' ;
* NOM         : PRPCB
*                Conductivité thermique d'un combustible CerCer
* Input : T ou TS : temperature sous forme de CHPOINT ou FLOTTANT
*                   (en K)
*         Domaine de validité : 500 a 1300 degres C
*         PORO : Porosité de taille contrôlée 
*         BETA : Coefficient pour la prise en compte de la porosité
*                dans la loi de Maxwell-Eucken
*         IRRAD : Combustible irradié ou pas ?
*                 Si irradié, on choisit un taux de combustion de 10 at%
*                 pour le (U,PU)C et une température d'irradiation
*                 égale à 1125°C pour le SiC .
*         FVOL : Fraction volumique d'inclusions (U,Pu)C 
*                Le reste est de la matrice (SiC)
* Output : LUPU : conductivité thermique du (U,Pu) de meme type que T
*          LSIC : conductivité thermique du Sic de meme type que T
*          LCB  : conductivité thermique du combustible CerCer
'DEBPROC' PRPCB PORO*'FLOTTANT' BETA*'FLOTTANT'
                IRRAD*'LOGIQUE' FVOL*'FLOTTANT'
                T/'CHPOINT' TS/'FLOTTANT' ;
* (U,Pu)C
C1 = 12.76D0 ;
C2 = 8.71D-3 ;
C3 = -1.88D-6 ;
* SiC
A = 42.58D0 ;
B = -1.5564D4 ;
C = 1.2977D7 ;
D = -1.8458D9 ;
* SiC irradié
E = 3.91112D-2 ;
F = 2.24732D-3 ;
TIRRADC = 1125.D0 ;
'SI' ('EXISTE' T) ;
   TC = '-' T 273.15D0 ;
   TC2 = '**' TC 2 ;
* (U,Pu)C non irradie et 100% dense
   LUPU1 = '+' ('+' ('*' TC C2) ('*' TC2 C3))
               C1 ;
* (U,Pu)C non irradie poreux
   CPOR = '/' ('-' 1.D0 PORO) ('+' 1.D0 ('*' BETA PORO)) ;
   LUPU2 = '*' LUPU1 CPOR ;
* (U,Pu)C irradie poreux
   'SI' IRRAD ;
      CIRRAD = ('-' 1.D0 ('*' 0.02D0 10.D0)) ;
      LUPU = '*' LUPU2 CIRRAD ;
   'SINON' ;
      LUPU = LUPU2 ;
   'FINSI' ;         
* SiC non irradie et 100% dense
   IT  = 'INVERSE' T ;
   IT2 = '**' IT 2 ;
   IT3 = '*' IT2 IT ;
   LSIC1 = '+' ('+' ('*' IT B)
                    ('+' ('*' IT2 C) ('*' IT3 D)))
               A ;
* SiC irradie               
   'SI' IRRAD ;
      CIRRAD = '*' ('EXP' ('*' TIRRADC F))
                   E ;
      LSIC = '*' LSIC1 CIRRAD ;
   'SINON' ;
      LSIC = LSIC1 ;                                  
   'FINSI' ;
* CerCer
   LMAT  = LSIC ;
   LINCL = LUPU ;
   LMAT2 = '*' LMAT 2.D0 ;
   LMI = '-' LMAT LINCL ;
   NUMER = '+' LINCL
               ('+' LMAT2
                    ('*' LMI ('*' -2.D0 FVOL))) ;
   DENOM = '+' LINCL
               ('+' LMAT2
                    ('*' LMI FVOL)) ;
   LCB = '*' LMAT ('*' NUMER ('INVERSE' DENOM)) ;
   'RESPROC' LUPU LSIC LCB ;
'FINSI' ; 
'SI' ('EXISTE' TS) ;
   TSC = '-' TS 273.15D0 ;
   TSC2 = '**' TSC 2 ;
* (U,Pu)C non irradie et 100% dense
   LUPU1 = '+' ('+' ('*' TSC C2) ('*' TSC2 C3))
               C1 ;
* (U,Pu)C non irradie poreux
   CPOR = '/' ('-' 1.D0 PORO) ('+' 1.D0 ('*' BETA PORO)) ;
   LUPU2 = '*' LUPU1 CPOR ;
* (U,Pu)C irradie poreux
   'SI' IRRAD ;
      CIRRAD = ('-' 1.D0 ('*' 0.02D0 10.D0)) ;
      LUPU = '*' LUPU2 CIRRAD ;
   'SINON' ;
      LUPU = LUPU2 ;
   'FINSI' ;         
* SiC non irradie et 100% dense
   ITS  = '/' 1.D0 TS ;
   ITS2 = '**' ITS 2 ;
   ITS3 = '*' ITS2 ITS ;
   LSIC1 = '+' ('+' ('*' ITS B)
                    ('+' ('*' ITS2 C) ('*' ITS3 D)))
               A ;
* SiC irradie               
   'SI' IRRAD ;
      CIRRAD = '*' ('EXP' ('*' TIRRADC F))
                   E ;
      LSIC = '*' LSIC1 CIRRAD ;
   'SINON' ;
      LSIC = LSIC1 ;                                  
   'FINSI' ;
* CerCer
   LMAT  = LSIC ;
   LINCL = LUPU ;
   LMAT2 = '*' LMAT 2.D0 ;
   LMI = '-' LMAT LINCL ;
   NUMER = '+' LINCL
               ('+' LMAT2
                    ('*' LMI ('*' -2.D0 FVOL))) ;
   DENOM = '+' LINCL
               ('+' LMAT2
                    ('*' LMI FVOL)) ;
   LCB = '*' LMAT ('*' NUMER ('/' 1.D0 DENOM)) ;
   'RESPROC' LUPU LSIC LCB ;
'FINSI' ; 
'FINPROC' ;
* NOM         : CALCUL
* DESCRIPTION : Boucle de calcul
*     bouche = VRAI si le canal d'helium central est bouche
*            = FAUX sinon 
'DEBPROC' CALCUL ;
'ARGUMENT' bouche*'LOGIQUE ' ;
*
* Puissance volumique dans le combustible
*
voltot = haut '*' pas '*' pas '*' t30 '*' 0.5D0 ;
volcb = 'MAXIMUM' ('RESULT' ('DOMA' $mtcb 'VOLUME')) ;
pmoy = '*' pvol ('/' voltot volcb) ;
'SI' ('NON' pvar) ;
   pcb = pmoy ;
'SINON' ;
   gamma = 1.5D0 ;
   raddeg = '/' 180.D0 PI ;
   zcb = 'COORDONNEE' 3 ('DOMA' $mtcb 'CENTRE') ;
   zmoy = '/' haut 2.D0 ;
   zcb = '+' zcb ('*' zmoy -1.D0) ;
   pmax = '*' pmoy ('/' gamma ('SIN' ('*' raddeg gamma))) ;
   fccos = '/' (2.D0 '*' gamma '*' raddeg) haut ;
   pcb1 = '*' ('COS' ('*' zcb fccos)) pmax ;
   chvol = 'DOMA' $mtcb 'VOLUME' ;
   puis1 = 'MAXIMUM' ('RESULT' ('*' chvol pmoy)) ;
   puis2 = 'MAXIMUM' ('RESULT' ('*' chvol pcb1)) ;
   'SI' ('NON' ('EGA' puis1 puis2 ('*' 0.01 ('ABS' puis2)))) ;
      'MESSAGE' ('CHAINE' 'puis1 = ' puis1 ) ;
      'MESSAGE' ('CHAINE' 'puis2 = ' puis2 ) ;
      'ERREUR' 'Les deux puissances devraient etre a peu pres egales' ;
   'FINSI' ;         
* Correction
   facor = '/' puis1 puis2 ;
   pcb2 = '*' pcb1 facor ;
   pcb = 'KCHT' $mtcb 'SCAL' 'CENTRE' pcb2 ;
'FINSI' ;   
*
* Débit d'hélium
*
qhe = qhenom ;
'MESSAGE' ('CHAINE' 'Débit=' qhe ' kg.m-2.s-1') ;      
*
* Initialisation des tables pour les calculs thermiques
* dans le combustible et dans l'hélium
*
rvcb = 'EQEX' 'ITMA' 1 'NITER' 1 'FIDT' 1000 ;
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcb 'OPER' 'LAPN' 'LCB' 'INCO' 'TCB' 
         ;
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcb 'OPER' 'FIMP' 'PCB' 'INCO' 'TCB' 
         ;
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtif 'OPER' 'ECHI' 'H' 'T0' 'INCO' 'TCB' 'TCB' 
         ;
rvhe = 'EQEX' 'ITMA' 1 'NITER' 1 'FIDT' 1000 ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL' 'CENTREE' 'NOCONS'
         'ZONE' $mthe 'OPER' 'KONV' 'RCPHE' 'UHE' 'LHE' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'LAPN' 'LHE' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'FIMP' 'PHE' 'INCO' 'THE' 
         ;
rvhe = 'EQEX' rvhe
         'CLIM' 'THE' 'TIMP' ('DOMA' $bahe 'MAILLAGE') theek ;
rvi = 'TABLE' 'INCO' ;
rvcb . 'INCO' = rvi ;
rvhe . 'INCO' = rvi ;
rvi . 'TCB' = 'KCHT'  $mtcb 'SCAL' 'SOMMET' theek ;
lupu lsic lcb = PRPCB por bet irrad fvo theek ;
rvi . 'LCB' = 'KCHT' $mtcb 'SCAL' 'CENTRE' lcb ;
rvi . 'PCB' = pcb ;
rvi . 'T0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H'  = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'THE' = 'KCHT' $mthe 'SCAL' 'SOMMET' theek ;
rvi . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
*
* Initialisation des propriétés physiques de l'hélium
*
rhe muhe lhe cphe = PRPHE prhe theek ;
uhe = '/' qhe rhe ;
rvi . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhe cphe) ;
rvi . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhe ;
rvi . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' (0.D0 0.D0 uhe) ;
rvi . 'QHE' = qhe ;
***************************************
*
* BOUCLE DE CALCUL
*
***************************************
'REPETER' bouc 10 ;
*   1) Calcul de la température moyenne dans l'hélium
   tmoyhe = 'NOEL' $mthe (rvi . 'THE') ;
      'MESSAGE' 'Maxi Temperature moyenne helium (degres C) : '
             ('-' ('MAXIMUM' tmoyhe) 273.15D0) ;
*   2) Transfert sur l'interface
   tmtif  = '*' mcf tmoyhe ;
   rvi . 'T0' = 'KCHT' $mtif 'SCAL' 'CENTRE' tmtif ;
*   3) Calcul du coefficient d'échange à l'interface
*      et des propriétés physiques du combustible
*      à partir des températures à l'instant précédent.
* Température de mélange
   tm = rvi . 'T0' ;
* Température à la paroi
   tp = 'NOEL' $mtif (rvi . 'TCB') ;
   t0 = '*' ('+' tm tp) 0.5D0 ;
*   t0 = rvi . 'T0' ;
   rho0 mu0 l0 cp0 = PRPHE prhe t0 ;
   pran0 = '*' ('*' mu0 cp0) ('INVERSE' l0) ;
   reyn0 = '*' ('*' qhe diam) ('INVERSE' mu0) ;
* Corrélation de Colburn
   nuss0 = 0.023 '*' ('**' reyn0 0.8D0)
                 '*' ('**' pran0 ('/' 1.D0 3.D0)) ;
   h0 = '/' ('*' nuss0 l0) diam ;
   'SI' (bouche) ;
     h0b = 'KCHT' $mtifb 'SCAL' 'CENTRE' h0 ;
     h0a = 'KCHT' $mtifa 'SCAL' 'CENTRE' 0.D0 ;
     rvi . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' (h0a 'ET' h0b) ;
   'SINON' ;
     rvi . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' h0 ;
   'FINSI' ;   
* Propriétés physiques du combustible
   tcbck = 'NOEL' $mtcb (rvi . 'TCB') ;
   lupu lsic lcb = PRPCB por bet irrad fvo tcbck ;
   rvi . 'LCB' = 'KCHT' $mtcb 'SCAL' 'CENTRE' lcb ;
*   4) Résolution de la thermique dans le combustible 
   EXEC rvcb ;
*   5) Calcul du flux de chaleur à la surface du combustible
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif 'OPER' 'MDIA' 'H' 'INCO' 'TCB' 'F0' ;
   rvmdia1 . 'INCO' = rvi ;                   
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif 'OPER' 'MDIA' 'H' 'INCO' 'T0'  'F0' ;
   rvmdia2 . 'INCO' = rvi ;                   
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f0 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T0'  (rvi . 'T0'))))
                ;
      'MESSAGE' 'Flux integre sur la paroi (W) : '
             ('MAXIMUM' ('RESULT' f0)) ;
*   6) Transfert en une puissance volumique dans l'hélium
   rvi . 'F0' = f0 ;
   phe = '*' mfc f0 ;
   rvi . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe ;
*
*   7) Calcul des propriétés de physiques de l'hélium 
*      a partir des propriétés à l'instant précédent.
*
   theck = 'NOEL' $mthe (rvi . 'THE') ;
   thesk = rvi . 'THE' ;
   rhec muhec lhec cphec = PRPHE prhe theck ;
   rhes muhes lhes cphes = PRPHE prhe thesk ;
   uhes = 'NOMC' 'UZ' ('*' qhe ('INVERSE' rhes)) ;
   rvi . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhec cphec) ;
   rvi . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhec ;
   rvi . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhes ;
*   8) Calcul de la température dans l'hélium
   EXEC rvhe ;
'FIN' bouc ;
'RESPRO' rvi ;
*
* End of procedure file CALCUL
*
'FINPROC' ;
*
* Fin des procédures
*
************************************************************************
'OPTION' 'ECHO' 1 ;
'TEMPS' 'ZERO' ;
'OPTION' 'DIME' 3 ;
'OPTION' 'ISOV' 'SULI' ;
nbisov = 25 ;
'SI' ('NON' interact) ;
  'OPTION' 'TRAC' 'PSC' ;
   'OPTION' 'ECHO' 0 ;
'SINON' ;
  'OPTION' 'TRAC' 'X' ;
'FINSI' ;  
formflot ='(1PE9.2)' ;
***************************************
*
* MAILLAGE
*
***************************************
* géométrie
pas  = 21.1D-3 ;
diam = 16.5D-3 ;
haut  = 2.D0 ;
*nbase = 10 ;
*nhaut = 20 ;
nbase = 7 ;
nhaut = 15 ;
*
* Physique
*
* gaz
prhe = 90.D0 ;
* débit d'hélium (kg.m-2.s-1)
qhenom = 210.D0 ;
** en K !
theec = 480.D0 ;
theek = '+' theec 273.15D0 ;
*
* Combustible
*
*
* Porosité (15%) et coefficient beta (2)
*
por = 0.15D0 ;
bet = 2.D0 ;
*
* Fraction volumique d'inclusions (U,Pu)C (70% pour
* le combustible performant)
*
fvo = 0.7D0 ;
*
* Attention a la définition de la puissance volumique.
* Ici, c'est la puissance totale dégagée par le combustible
* divisée par le volume total (combustible+refroidissement)
* Donc, il y a une correction pour avoir la puissance générée par unité
* de volume de combustible (cf. variable pvol)
*
pvol = 1.D+8 ;
titglob = 'CHAINE' ' ' ;
'SI' irrad ;
   titglob = 'CHAINE' titglob 'Cb irradie' ;
'SINON' ;   
   titglob = 'CHAINE' titglob 'Cb neuf' ;
'FINSI' ;   
q = '/' pas 2.D0 ;
R = '/' diam 2.D0 ;
s30 = 'SIN' 30.D0 ;
c30 = 'COS' 30.D0 ;
t30 = '/' s30 c30 ;
* longueurs de chaque ligne
lcera = ('*' ('/' PI 6.D0) R) ;
lcerb = ('*' ('/' PI 2.D0) R) ;
lpdia = '-' pas diam ;
lgdia = '-' ('/' pas c30) R ;
lhex  = '-' ('*' pas t30) R ;
* On choisit comme échelle la plus petite des longueurs
lbase = 'MINIMUM' ('PROG' lcera lcerb lpdia lgdia lhex) ;
* On calcule les facteurs qui permettent d'avoir un maillage
* a peu près régulier en longueurs
fcera = '/' lcera lbase ;
fcerb = '/' lcerb lbase ;
fpdia = '/' lpdia lbase ;
fgdia = '/' lgdia lbase ;
fhex  = '/' lhex lbase ;
* maillage
rbase = 'FLOTTANT' nbase ;
ncera  = '*' rbase fcera ;
ncerb  = '*' rbase fcerb ;
nhex  = '*' rbase fhex ;
npdia  = '*' rbase fpdia ;
ngdia  = '*' rbase fgdia ;
*
ncera = 'ENTIER' ('+' ncera 0.99D0) ;
ncerb = 'ENTIER' ('+' ncerb 0.99D0) ;
nhex = 'ENTIER' ('+' nhex 0.99D0) ;
npdia = 'ENTIER' ('+' npdia 0.99D0) ;
ngdia = 'ENTIER' ('+' ngdia 0.99D0) ;
'SI' ('EGA' discr 'LINE') ;
   'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
'SINON' ;   
   'OPTION' 'DIME' 2 'ELEM' 'QUA8' ;
'FINSI' ;
* définition des points pour le combustible
p0 = 0.D0 0.D0 ;
p1 = R    0.D0 ;
p2 = ('-' pas R) 0.D0 ;
p3 = pas  0.D0 ;
p4 = pas  R    ;
p5 = pas  ('*' pas t30) ;
p6 = ('*' R c30) ('*' R s30) ;
* définition des lignes pour le combustible
l1 = 'DROIT' npdia p1 p2 ;
l2 = 'CERCLE' ncerb p2 p3 p4 ;
l3 = 'DROIT' nhex  p4 p5 ;
l4 = 'DROIT' ngdia p5 p6 ;
l5 = 'CERCLE' ncera p6 p0 p1 ;
* définition des surfaces pour le combustible
ctn  = l1 'ET' l2 'ET' l3 'ET' l4 'ET' l5 ;
stcb = 'SURFACE' ctn 'PLANE' ;
* définition des points pour l'hélium (deux  triangles isocèles de
* surfaces égales aux secteurs circulaires)
Rpa = '*' R ('**' ('/' PI 3.D0) 0.5D0) ;
p7 = Rpa   0.D0 ;
p8 = ('*' Rpa c30) ('*' Rpa s30) ;
sthea = 'MANUEL' 'TRI3' p0 p7 p8 ;
Rpb = '*' R ('**' ('/' PI 2.D0) 0.5D0) ;
p9  = ('-' pas Rpb) 0.D0 ;
p10 = pas Rpb ;
stheb = 'MANUEL' 'TRI3' p9 p3 p10 ;
'SI' ('NON' ('EGA' discr 'LINE')) ;
   sthea = 'CHANGER' 'QUADRATIQUE' sthea ;
   stheb = 'CHANGER' 'QUADRATIQUE' stheb ;
'FINSI' ;
*'TRACER' (stcb 'ET' sthea 'ET' stheb) ;
*
'SI' ('EGA' discr 'LINE') ;
   'OPTION' 'DIME' 3 'ELEM' 'CUB8' ;
'SINON' ;   
   'OPTION' 'DIME' 3 'ELEM' 'CU20' ;
'FINSI' ;
hnh = '/' haut ('FLOTTANT' nhaut) ;
vechau = 0.D0 0.D0 hnh ;
* définition des éléments dans une hauteur élémentaire
bahea = sthea ;
baheb = stheb ;
mthea1 = 'VOLUME' sthea 1 'TRAN' vechau ;
mtheb1 = 'VOLUME' stheb 1 'TRAN' vechau ;
mtifa1 = 'TRANSLATION' l5   1   vechau ;
mtifb1 = 'TRANSLATION' l2   1   vechau ;
mtcb1 = 'VOLUME' stcb 1 'TRAN' vechau ;
   mthea = mthea1 ;
   mtheb = mtheb1 ;
   mtifa = mtifa1 ;
   mtifb = mtifb1 ;
   mtcb = mtcb1 ;
tabfac = 'TABLE' ;
tabele = 'TABLE' ;
icpt = 1 ;
nmtifa1 = 'NBEL' mtifa1 ;
'REPETER' imtifa1 nmtifa1 ;
   lmtifa1 = 'ELEM' mtifa1 &imtifa1 ;
   tabfac . icpt = lmtifa1 ;
   tabele . icpt = mthea1 ;
   icpt = '+' icpt 1 ;
'FIN' imtifa1 ;
nmtifb1 = 'NBEL' mtifb1 ;
'REPETER' imtifb1 nmtifb1 ;
   lmtifb1 = 'ELEM' mtifb1 &imtifb1 ;
   tabfac . icpt = lmtifb1 ;
   tabele . icpt = mtheb1 ;
   icpt = '+' icpt 1 ;
'FIN' imtifb1 ;
* Génération par translation de la hauteur
'REPETER' ihaut ('-' nhaut 1) ;
   vectran = 0.D0 0.D0 ('*' hnh ('FLOTTANT' &ihaut)) ;
   mthea2 mtifa2 mtheb2 mtifb2 mtcb2 =
      'PLUS' mthea1 mtifa1 mtheb1 mtifb1 mtcb1 vectran ;
   nmtifa2 = 'NBEL' mtifa2 ;
   'REPETER' imtifa2 nmtifa2 ;
      lmtifa2 = 'ELEM' mtifa2 &imtifa2 ;
      tabfac . icpt = lmtifa2 ;
      tabele . icpt = mthea2 ;
      icpt = '+' icpt 1 ;
   'FIN' imtifa2 ;
   nmtifb2 = 'NBEL' mtifb2 ;
   'REPETER' imtifb2 nmtifb2 ;
      lmtifb2 = 'ELEM' mtifb2 &imtifb2 ;
      tabfac . icpt = lmtifb2 ;
      tabele . icpt = mtheb2 ;
      icpt = '+' icpt 1 ;
   'FIN' imtifb2 ;
   mthea = 'ET' mthea2 mthea ;
   mtifa = 'ET' mtifa2 mtifa ;
   mtheb = 'ET' mtheb2 mtheb ;
   mtifb = 'ET' mtifb2 mtifb ;
   mtcb = 'ET' mtcb2 mtcb ;
'FIN' ihaut ;   
bahe = 'ET' bahea baheb ;
mthe = 'ET' mthea mtheb ;
mtif = 'ET' mtifa mtifb ;
* Génération des lignes dans l'hélium
vechaut = 0.D0 0.D0 haut ;
lighea1 = p0 'DROIT' nhaut ('PLUS' p0 vechaut) ;
lighea2 = p7 'DROIT' nhaut ('PLUS' p7 vechaut) ;
lighea3 = p8 'DROIT' nhaut ('PLUS' p8 vechaut) ;
ligheb1 = p9 'DROIT' nhaut ('PLUS' p9 vechaut) ;
ligheb2 = p3 'DROIT' nhaut ('PLUS' p3 vechaut) ;
ligheb3 = p10 'DROIT' nhaut ('PLUS' p10 vechaut) ;
lighea  = lighea1 'ET' lighea2 'ET' lighea3 ;
ligheb  = ligheb1 'ET' ligheb2 'ET' ligheb3 ;
lighe = 'ET' lighea ligheb ;
mt = ('ET' mthe mtcb) ;
pelim = '*' ('**' ('/' ('MESURE' mt) ('NBEL' mt))
                  ('/' 1.D0 ('VALEUR' 'DIME')))
            1.D-3 ;
'ELIMINATION' (mtif 'ET' mtcb) pelim ;
'ELIMINATION' (mthe 'ET' bahe 'ET' lighe)
              pelim ;
*
* Génération des QUAFs
*
_bahe = 'CHANGER' bahe 'QUAF' ;
_bahea = 'CHANGER' bahea 'QUAF' ;
_baheb = 'CHANGER' baheb 'QUAF' ;
_mthe = 'CHANGER' mthe 'QUAF' ;
_mthea = 'CHANGER' mthea 'QUAF' ;
_mtheb = 'CHANGER' mtheb 'QUAF' ;
_mtif = 'CHANGER' mtif 'QUAF' ;
_mtifa = 'CHANGER' mtifa 'QUAF' ;
_mtifb = 'CHANGER' mtifb 'QUAF' ;
_mtcb = 'CHANGER' mtcb 'QUAF' ;
***************************************
*
* CREATION DE LA MATRICE DE COUPLAGE
*
***************************************
*
* Génération de connectivites face-centre cfc 
* et des matrices de passage centre-face pour projeter T0
* et face-centre pour ramener une intégrale de flux
* a une puissance volumique
*
nitab = 'DIME' tabfac ;
'REPETER' iitab nitab ;
   ifac = tabfac . &iitab ;
   iqfac = 'CHANGER' ifac 'QUAF' ;
   miqfac = 'MODELISER' iqfac 'NAVIER_STOKES' 'QUAF';
   ciqfac = ('DOMA' miqfac 'CENTRE') 'POIN' 1 ;
   iele = tabele . &iitab ;
   iqele = 'CHANGER' iele 'QUAF' ;
   miqele = 'MODELISER' iqele 'NAVIER_STOKES' 'QUAF';
   ciqele = ('DOMA' miqele 'CENTRE') 'POIN' 1 ;
* Attention ('MESU' iele) ne fonctionne pas en axisymetrique   
   volele = 'MAXIMUM' ('DOMA' miqele 'VOLUME') ;
   clfc = 'MANUEL' 'SEG2' ciqfac ciqele ;
   mlcf = 'MANUEL' 'RIGIDITE' clfc ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 1.D0) ('PROG' 0.D0 0.D0) ;
   mlfc = 'MANUEL' 'RIGIDITE' clfc ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 0.D0)
            ('PROG' ('/' 1.D0 volele) 0.D0) ;
   'SI' ('EGA' &iitab 1) ;
      cfc = clfc ;
      mcf = mlcf ;
      mfc = mlfc ;
   'SINON' ;
      cfc = 'ET' cfc clfc ;
      mcf = 'ET' mcf mlcf ;
      mfc = 'ET' mfc mlfc ;
   'FINSI' ;
'FIN' iitab ;
'ELIMINATION' (_bahe 'ET' _bahea 'ET' _baheb
          'ET' _mthe 'ET' _mthea 'ET' _mtheb 'ET' cfc) pelim ;
'ELIMINATION' (_mtcb 'ET'
               _mtif 'ET' _mtifa 'ET' _mtifb 'ET' cfc) pelim ;
mcf = 'CHANGER' mcf 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc = 'CHANGER' mfc 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
***************************************
*
* PREPARATION DU CALCUL
*
***************************************
*
* Modèles Navier-Stokes
*
$bahe = 'MODELISER' _bahe 'NAVIER_STOKES' discr ;
$bahea = 'MODELISER' _bahea 'NAVIER_STOKES' discr ;
$baheb = 'MODELISER' _baheb 'NAVIER_STOKES' discr ;
$mthe = 'MODELISER' _mthe 'NAVIER_STOKES' discr ;
$mthea = 'MODELISER' _mthea 'NAVIER_STOKES' discr ;
$mtheb = 'MODELISER' _mtheb 'NAVIER_STOKES' discr ;
$mtif = 'MODELISER' _mtif 'NAVIER_STOKES' discr ;
$mtifa = 'MODELISER' _mtifa 'NAVIER_STOKES' discr ;
$mtifb = 'MODELISER' _mtifb 'NAVIER_STOKES' discr ;
$mtcb = 'MODELISER' _mtcb 'NAVIER_STOKES' discr ;
rvi = CALCUL FAUX ;
rwi = CALCUL VRAI ;
tcpu = 'TEMPS' 'NOEC' ;
tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
'MESSAGE' ('CHAINE' 'tcpus=' tcpus) ;
***************************************
*
* POST-TRAITEMENT
*
***************************************
'SI' graph ;
*
* Maillage (+ changement d'échelle adéquat)
* 
   _mt = 'ET' _mthe _mtcb ;
   mt = 'ET' mthe mtcb ;
   zmt = 'COORDONNEE' 3 _mt ;
   factech =  ('-' ('/' 1.D0 75.D0) 1.D0) ;
   dxmt = 'NOMC' 'UZ' ('*' zmt factech) 'NATURE' 'DISCRET' ;
   orig = 'FORME' ;
   'FORME' dxmt ;
   xmthea ymthea zmthea = 'COORDONNEE' mthea ;
   factech = ('-' ('/' 1.D0 4.D0) 1.D0) ;
   dxmta = 'ET'
             ('NOMC' 'UX' ('*' xmthea factech) 'NATURE' 'DISCRET')
             ('NOMC' 'UY' ('*' ymthea factech) 'NATURE' 'DISCRET')
             ;
   xmtheb ymtheb zmtheb = 'COORDONNEE' mtheb ;
   dxmtb = 'ET'
             ('NOMC' 'UX' ('*' ('-' xmtheb pas) factech)
              'NATURE' 'DISCRET')
             ('NOMC' 'UY' ('*' ymtheb factech) 'NATURE' 'DISCRET')
             ;
   dxmt = dxmta 'ET' dxmtb ;          
* manipulation maillage et champ de temperature
   tcbn = ('-' (rvi . 'TCB') 273.15D0) ;
   tcbb = ('-' (rwi . 'TCB') 273.15D0) ;
   tcbn2 tcbb2 mtcb2 = tcbn tcbb mtcb 'TOURNER' 150.D0
      (0.D0 0.D0 0.D0) (0.D0 0.D0 1.D0) ;
   then = ('-' (rvi . 'THE') 273.15D0) ;
   theb = ('-' (rwi . 'THE') 273.15D0) ;
   then2 theb2 mthe2 = then theb mthe 'PLUS' dxmt ;
   xmt ymt zmt = 'COORDONNEE' _mt ;
   mxmt = 'MAXIMUM' xmt ;
   mymt = 'MAXIMUM' ymt ;
   mzmt = 'MAXIMUM' zmt ;
   oeil = (('*' mxmt 3.D0) ('*' mymt -4.D0)  ('*' mzmt 2.85D0)) ; 
   oeil2 = (('*' mxmt 3.D0) ('*' mymt 4.D0)  ('*' mzmt 1.85D0)) ; 
   mail = mtcb 'ET' mthe2 ;
   npomail = 'NBNO' mail ;
   nelmail = 'NBEL' mail ;
   titmail = 'CHAINE' 'Maillage ' 'NBPO=' npomail
                                ' NBELEM=' nelmail titglob ;
   'TRACER' oeil 'CACH' mail 'TITR' titmail ;
*               
* Tracé des températures, du coefficient d'échange et de la
* composante verticale des vitesses
*
   'REPETER' bctem 3 ;
      'SI' ('EGA' &bctem 1) ;
         tcbc  = tcbn ;
         tcbc2 = tcbn2 ;
         thec  = then ;
         thec2 = then2 ;
         h  = rvi . 'H' ;
         uz = 'EXCO' 'UZ' (rvi . 'UHE') ;
         tit2  = 'CHAINE' ' Nominal' ;
      'FINSI' ;   
      'SI' ('EGA' &bctem 2) ;
         tcbc  = tcbb ;
         tcbc2 = tcbb2 ;
         thec  = theb ;
         thec2 = theb2 ;
         h  = rwi . 'H' ;
         uz = 'EXCO' 'UZ' (rwi . 'UHE') ;
         tit2  = 'CHAINE' ' Accidentel' ;
      'FINSI' ;   
      'SI' ('EGA' &bctem 3) ;
         tcbc  = '-' tcbb tcbn ;
         tcbc2 = '-' tcbb2 tcbn2 ;
         thec  = '-' theb then ;
         thec2 = '-' theb2 then2 ;
         h  = '-' (rwi . 'H') (rvi . 'H') ;
         uz = '-' ('EXCO' 'UZ' (rwi . 'UHE'))
                  ('EXCO' 'UZ' (rvi . 'UHE')) ;
         tit2  = 'CHAINE' ' Acc. - Nom.' ;
      'FINSI' ;   
      tit = 'CHAINE' 'Temp. (°C)' tit2 ;
      tracchpo oeil (tcbc 'ET' tcbc2 'ET' thec2)
                    (mtcb 'ET' mtcb2 'ET' mthe2)
                    tit ;
      tit = 'CHAINE' 'Temp. Comb. (°C)' tit2 ;
*   tracchpo oeil tcbc mtcb tit ;
*   tracchpo oeil tcbc2 mtcb2 tit ;
      tracchpo oeil (tcbc 'ET' tcbc2) (mtcb 'ET' mtcb2) tit ;
*   tracchpo oeil tcbc2 mtcb2 tit ;
* Coefficient d'echange
      tit = 'CHAINE' 'Coeff. echange (W.m-2.K-1)' tit2 ;
      tracchml oeil $mtif h mtif tit ;
      tit = 'CHAINE' 'Temp. He (°C)' tit2 ;
      tracchpo oeil thec mthe tit ;
      tit = 'CHAINE' 'Uz He (m.s-1)' tit2 ;
      tracchpo oeil uz  mthe tit ;
   'FIN' bctem ;
   'FORME' orig ;
*
* Evolutions temperature et vitesse
* dans les deux tubes d'hélium nominal et bouche
*
   evthean = 'EVOL' 'CHPO' then lighea1 ;
   evthebn = 'EVOL' 'CHPO' then ligheb1 ;
   evtheab = 'EVOL' 'JAUN' 'CHPO' theb lighea1 ;
   evthebb = 'EVOL' 'JAUN' 'CHPO' theb ligheb1 ;
   evthe = evthean 'ET' evthebn 'ET' evtheab 'ET' evthebb ;
   titthe = 'CHAINE' 'Température Helium' titglob ;
   titxthe = 'CHAINE' 's (m)' ;
   titythe = 'CHAINE' 'T (°C)' ;
   tthe = 'TABLE' ;
   tthe . 1 = 'CHAINE' 'MARQ CROI' ;
   tthe . 2 = 'CHAINE' 'MARQ PLUS' ;
   tthe . 3 = 'CHAINE' 'MARQ CROI' ;
   tthe . 4 = 'CHAINE' 'MARQ PLUS' ;
   tthe . 'TITRE' = 'TABLE' ;
   tthe . 'TITRE' . 1 = 'CHAINE' 'Central nominal' ;
   tthe . 'TITRE' . 2 = 'CHAINE' 'Externe nominal' ;
   tthe . 'TITRE' . 3 = 'CHAINE' 'Central accidentel' ;
   tthe . 'TITRE' . 4 = 'CHAINE' 'Externe accidentel' ;
   'DESSIN' evthe 'TITR' titthe 'TITX' titxthe 'TITY' titythe
                'LEGE' 'MIMA' tthe ;
*
   uhen = 'EXCO' 'UZ' (rvi . 'UHE') ;
   uheb = 'EXCO' 'UZ' (rwi . 'UHE') ;
   evuhean = 'EVOL' 'CHPO' uhen lighea1 ;
   evuhebn = 'EVOL' 'CHPO' uhen ligheb1 ;
   evuheab = 'EVOL' 'JAUN' 'CHPO' uheb lighea1 ;
   evuhebb = 'EVOL' 'JAUN' 'CHPO' uheb ligheb1 ;
   evuhe = evuhean 'ET' evuhebn 'ET' evuheab 'ET' evuhebb ;
   tituhe = 'CHAINE' 'Vitesse verticale Helium' titglob ;
   titxuhe = 'CHAINE' 's (m)' ;
   tityuhe = 'CHAINE' 'uz (m.s-1)' ;
   tuhe = 'TABLE' ;
   tuhe . 1 = 'CHAINE' 'MARQ CROI' ;
   tuhe . 2 = 'CHAINE' 'MARQ PLUS' ;
   tuhe . 3 = 'CHAINE' 'MARQ CROI' ;
   tuhe . 4 = 'CHAINE' 'MARQ PLUS' ;
   tuhe . 'TITRE' = 'TABLE' ;
   tuhe . 'TITRE' . 1 = 'CHAINE' 'Central nominal' ;
   tuhe . 'TITRE' . 2 = 'CHAINE' 'Externe nominal' ;
   tuhe . 'TITRE' . 3 = 'CHAINE' 'Central accidentel' ;
   tuhe . 'TITRE' . 4 = 'CHAINE' 'Externe accidentel' ;
   'DESSIN' evuhe 'TITR' tituhe 'TITX' titxuhe 'TITY' tityuhe
                   'LEGE' 'MIMA' tuhe ;
'FINSI' ;                   
*
* Tableau recapitulatif et vérifications diverses
*
optecho = 'VALEUR' 'ECHO' ;
testerr = FAUX ;
'OPTION' 'ECHO' 0 ;
'MESSAGE' '**************' ;
'MESSAGE' ('CHAINE' 'couplage_TH1D_Th3D_2 :' titglob) ;
mail = mtcb 'ET' mthe ;
npomail = 'NBNO' mail ;
nelmail = 'NBEL' mail ;
'MESSAGE' ('CHAINE' 'nb. points=' npomail ' ; nb. elem=' nelmail) ;
'MESSAGE' ('CHAINE' 'Temps CPU (s)=' tcpus) ;
'MESSAGE' '-----------' ;
'MESSAGE' 'Verifications :' ;
'REPETER' bclsit 2 ;
   'SI' ('EGA' &bclsit 1) ;
      thec = rvi . 'THE' ;
      uhec = rvi . 'UHE' ;
      titsit = 'CHAINE' 'Nominal' ;
   'FINSI' ;
   'SI' ('EGA' &bclsit 2) ;
      thec = rwi . 'THE' ;
      uhec = rwi . 'UHE' ;
      titsit = 'CHAINE' 'Accidentel' ;
   'FINSI' ;
      'MESSAGE' ('CHAINE' '--> ' titsit) ;
   'REPETER' bcltub 2 ;
      'SI' ('EGA' &bcltub 1) ;
         lighe1 = lighea1 ;
         lighe2 = lighea2 ;
         lighe3 = lighea3 ;
         tittub = 'CHAINE' 'Canal Central' ;
      'FINSI' ;
      'SI' ('EGA' &bcltub 2) ;
         lighe1 = ligheb1 ;
         lighe2 = ligheb2 ;
         lighe3 = ligheb3 ;
         tittub = 'CHAINE' 'Canal Externe' ;
      'FINSI' ;
         'MESSAGE' ('CHAINE' '----> ' tittub) ;
* Symétrie température dans lhélium
      evthe1 = 'EVOL' 'CHPO' thec lighe1 ;
      evthe2 = 'EVOL' 'CHPO' thec lighe2 ;
      evthe3 = 'EVOL' 'CHPO' thec lighe3 ;
      the1 = 'EXTRAIRE' evthe1 'ORDO' ;
      the2 = 'EXTRAIRE' evthe2 'ORDO' ;
      the3 = 'EXTRAIRE' evthe3 'ORDO' ;
      the2m1 = 'MAXIMUM' ('-' the2 the1) 'ABS' ;
      the3m1 = 'MAXIMUM' ('-' the3 the1) 'ABS' ;
      testsym1 = > the2m1 1.D-5 ;
      testsym2 = > the3m1 1.D-5 ;
         'MESSAGE' 'Symétrie température dans lhélium :' ;
         'MESSAGE' ('CHAINE' 'the2m1=' the2m1) ;
         'MESSAGE' ('CHAINE' 'the3m1=' the3m1) ;
      'SI' ('OU' testsym1 testsym2) ;
         'MESSAGE' '!!Temperature non symetrique dans lhelium' ;
      'FINSI' ;
      testerr = testerr 'OU' testsym1 'OU' testsym2 ;
* Symétrie Uz dans lhélium
      evuzhe1 = 'EVOL' 'CHPO' uhec 'UZ' lighe1 ;
      evuzhe2 = 'EVOL' 'CHPO' uhec 'UZ' lighe2 ;
      evuzhe3 = 'EVOL' 'CHPO' uhec 'UZ' lighe3 ;
      uzhe1 = 'EXTRAIRE' evuzhe1 'ORDO' ;
      uzhe2 = 'EXTRAIRE' evuzhe2 'ORDO' ;
      uzhe3 = 'EXTRAIRE' evuzhe3 'ORDO' ;
      uzhe2m1 = 'MAXIMUM' ('-' uzhe2 uzhe1) 'ABS' ;
      uzhe3m1 = 'MAXIMUM' ('-' uzhe3 uzhe1) 'ABS' ;
      testsym1 = > uzhe2m1 1.D-5 ;
      testsym2 = > uzhe3m1 1.D-5 ;
         'MESSAGE' 'Symétrie Uz dans lhélium :' ;
         'MESSAGE' ('CHAINE' 'uzhe2m1=' uzhe2m1) ;
         'MESSAGE' ('CHAINE' 'uzhe3m1=' uzhe3m1) ;
      'SI' ('OU' testsym1 testsym2) ;
         'MESSAGE' '!!Uz non symetrique dans lhelium' ;
      'FINSI' ;
      testerr = testerr 'OU' testsym1 'OU' testsym2 ;
   'FIN' bcltub ;
* Nullité Vitesse horizontale dans lhélium
   uhex = 'EXCO' uhec 'UX' ;
   uhey = 'EXCO' uhec 'UY' ;
   muhex = 'MAXIMUM' uhex 'ABS' ;
   muhey = 'MAXIMUM' uhey 'ABS' ;
   testsymx = > muhex 1.D-5 ;
   testsymy = > muhey 1.D-5 ;
      'MESSAGE' 'Nullité Vitesse horizontale dans lhelium :' ;
      'MESSAGE' ('CHAINE' 'maxuhex=' muhex) ;
      'MESSAGE' ('CHAINE' 'maxuhey=' muhey) ;
   'SI' ('OU' testsymx testsymy) ;
      'MESSAGE' '!!Composante UX ou UY non nulle dans lhelium' ;
   'FINSI' ;
   testerr = testerr 'OU' testsymx 'OU' testsymy ;
'FIN' bclsit ;   
*
* Verification des bilans d'énergie
*
'MESSAGE' '-----------' ;
'MESSAGE' 'Bilans denergie :' ;
'REPETER' bclsit 2 ;
   'SI' ('EGA' &bclsit 1) ;
      rti = rvi ;
      titsit = 'CHAINE' 'Nominal' ;
   'FINSI' ;
   'SI' ('EGA' &bclsit 2) ;
      rti = rwi ;
      titsit = 'CHAINE' 'Accidentel' ;
   'FINSI' ;
   'MESSAGE' ('CHAINE' '--> ' titsit) ;
*  Puissance dégagée dans le combustible
   pcb = rti . 'PCB' ;
   chvolcb = 'DOMA' $mtcb 'VOLUME' ;
   pdegcb = 'MAXIMUM' ('RESULT' ('*' chvolcb pcb)) ;
*  Flux intégré sur la paroi
   f0 = rti . 'F0' ;
   pdegif = 'MAXIMUM' ('RESULT' f0) ;
   errpif = '*' (ERRREL pdegif pdegcb) 100.D0 ;
*  Puissance acquise par l'hélium
*   Température en sortie (en K)
   the = rti . 'THE' ;
   evthea = 'EVOL' 'CHPO' the lighea1 ;
   thea = 'EXTRAIRE' evthea 'ORDO' ;
   tsora = 'EXTRAIRE' thea ('DIME' thea) ;
   tenta = 'EXTRAIRE' thea 1 ;
   shea  = 'MAXIMUM' ('RESULT' ('DOMA' $bahea 'VOLUME')) ;
   evtheb = 'EVOL' 'CHPO' the ligheb1 ;
   theb = 'EXTRAIRE' evtheb 'ORDO' ;
   tsorb = 'EXTRAIRE' theb ('DIME' theb) ;
   tentb = 'EXTRAIRE' theb 1 ;
   sheb  = 'MAXIMUM' ('RESULT' ('DOMA' $baheb 'VOLUME')) ;
*   On utilise le fait que le cp de l'hélium
*   est constant
   rhe muhe lhe cphe = PRPHE prhe tenta ;
   qhe = rti . 'QHE' ;
   pdeghea = qhe '*' shea '*' cphe '*' ('-' tsora tenta) ;
   pdegheb = qhe '*' sheb '*' cphe '*' ('-' tsorb tentb) ;
   pdeghe = '+' pdeghea pdegheb ;
   errphe = '*' (ERRREL pdeghe pdegcb) 100.D0 ;
   'MESSAGE' ('CHAINE' 'pdegcb=' pdegcb) ;
   'MESSAGE' ('CHAINE' 'pdegif=' pdegif
                       ' ; errpif=' errpif ' %') ;
   'MESSAGE' ('CHAINE' 'pdeghe=' pdeghe
                       ' ; errphe=' errphe ' %') ;
* Ce sont des erreurs en % (i.e. limite à 0.1%)
   testene1 = > errpif 1.D-1 ;
   testene2 = > errphe 1.D-1 ;
   'SI' ('OU' testene1 testene2) ;
      'MESSAGE' '!!Pb. conservation de lenergie' ;
   'FINSI' ;
   testerr = testerr 'OU' testene1 'OU' testene2 ;
* Verification du péclet de maille dans l'helium
   vmoyhe = 'NOEL' $mthe ('EXCO' 'UZ' (rti . 'UHE')) ;
   lmoyhea = '/' ('DOMA' $mthea 'VOLUME') shea ;
   lmoyheb = '/' ('DOMA' $mtheb 'VOLUME') sheb ;
   lmoyhe = lmoyhea 'ET' lmoyheb ;
   alphhe = '/' (rti. 'LHE') (rti . 'RCPHE') ;
   peclhe = '/' ('*' vmoyhe lmoyhe) alphhe ;
   pemin = 'MINIMUM' peclhe ;
   pemax = 'MAXIMUM' peclhe ;
   'MESSAGE' 'Peclet de maille helium :' ;
   'MESSAGE' ('CHAINE' 'pemin=' pemin) ;
   'MESSAGE' ('CHAINE' 'pemax=' pemax) ;
'FIN' bclsit ;   
* Temperature de sortie helium et temp. max combustible
'MESSAGE' '-----------' ;
'MESSAGE' 'Resultats Helium :' ;
'REPETER' bclsit 2 ;
   'SI' ('EGA' &bclsit 1) ;
      rti = rvi ;
      titsit = 'CHAINE' 'Nominal' ;
   'FINSI' ;
   'SI' ('EGA' &bclsit 2) ;
      rti = rwi ;
      titsit = 'CHAINE' 'Accidentel' ;
   'FINSI' ;
   'MESSAGE' ('CHAINE' '--> ' titsit) ;
   the = rti . 'THE' ;
   uhe = 'EXCO' 'UZ' (rti . 'UHE') ;
   'REPETER' bcltub 2 ;
      'SI' ('EGA' &bcltub 1) ;
         lighe = lighea1 ;
         tittub = 'CHAINE' 'Canal Central' ;
      'FINSI' ;
      'SI' ('EGA' &bcltub 2) ;
         lighe = ligheb1 ;
         tittub = 'CHAINE' 'Canal Externe' ;
      'FINSI' ;
      'MESSAGE' ('CHAINE' '----> ' tittub) ;
      evthe = 'EVOL' 'CHPO' the lighe ;
      evuhe = 'EVOL' 'CHPO' uhe lighe ;
      othe = 'EXTRAIRE' evthe 'ORDO' ;
      ouhe = 'EXTRAIRE' evuhe 'ORDO' ;
      thee = 'EXTRAIRE' othe 1 ;
      thes = 'EXTRAIRE' othe ('DIME' othe) ;
      uhee = 'EXTRAIRE' ouhe 1 ;
      uhes = 'EXTRAIRE' ouhe ('DIME' ouhe) ;
      'MESSAGE' ('CHAINE' 'Vit.  entree helium  =' uhee ' m.s-1') ;
      'MESSAGE' ('CHAINE' 'Vit.  sortie helium  =' uhes ' m.s-1') ;
      'MESSAGE' ('CHAINE' 'Temp. entree helium  =' ('-' thee 273.15D0)
                 ' degres C') ;
      'MESSAGE' ('CHAINE' 'Temp. sortie helium  =' ('-' thes 273.15D0)
                 ' degres C') ;
   'FIN' bcltub ;
'FIN' bclsit ;
* Tests sur Tmax (02/12/03)
'MESSAGE' '-----------' ;
'MESSAGE' 'Resultats Combustible :' ;
'REPETER' bclsit 2 ;
   'SI' ('EGA' &bclsit 1) ;
      rti = rvi ;
      titsit = 'CHAINE' 'Nominal' ;
      'SI' irrad ;
         tmref = 1.15428D3 ; altref = 1.2D0 ;
      'SINON' ;
         tmref = 1.10115D3 ; altref = 1.2D0 ;
      'FINSI' ;
   'FINSI' ;
   'SI' ('EGA' &bclsit 2) ;
      rti = rwi ;
      titsit = 'CHAINE' 'Accidentel' ;
      'SI' irrad ;
         tmref = 1.59486D3 ; altref = 1.2D0 ;
      'SINON' ;
         tmref = 1.46246D3 ; altref = 1.2D0 ;
      'FINSI' ;
   'FINSI' ;
   'MESSAGE' ('CHAINE' '--> ' titsit) ;
   tcb = rti . 'TCB' ;                 
   tmincbc = '-' ('MINIMUM' tcb) 273.15D0 ;
   amincbc = 'MAXIMUM' ('COORDONNEE' 3 (tcb 'POIN' 'MINI' 'ABS')) ;
   tmaxcbc = '-' ('MAXIMUM' tcb) 273.15D0 ;
   amaxcbc = 'MAXIMUM' ('COORDONNEE' 3 (tcb 'POIN' 'MAXI' 'ABS')) ;
   'MESSAGE' ('CHAINE' 'Temp. min combustible=' tmincbc ' degres C') ;
   'MESSAGE' ('CHAINE' '          altitude   =' amincbc ' m') ;
   'MESSAGE' ('CHAINE' 'Temp. max combustible=' tmaxcbc ' degres C') ;
   'MESSAGE' ('CHAINE' '          altitude   =' amaxcbc ' m') ;
* Tests sur Tmax (02/12/03)
* On utilise une discrétisation linéaire et une puissance variable
   'SI' ('ET' ('EGA' discr 'LINE') pvar) ;
      errtm = '*' (ERRREL tmaxcbc tmref) 100.D0 ;
      erral = '*' (ERRREL amaxcbc altref) 100.D0 ;
      'MESSAGE' ('CHAINE' 'Temp. max cb ref     =' tmref ' degres C') ;
      'MESSAGE' ('CHAINE' '        alt. ref     =' altref ' m') ;
      testtmax = 'OU' (> errtm 1.D-1) (> erral 1.D-1) ;
      'SI' testtmax ;
         'MESSAGE' '!! Pb. sur Tmax' ;
      'FINSI' ;
      testerr = testerr 'OU' testtmax ;
   'FINSI' ;      
'FIN' bclsit ;   
'MESSAGE' '**************' ;
'OPTION' 'ECHO' optecho ;
*!
*! Error checking
*!
'SI' testerr ;
   'MESSAGE' 'Il y a eu des erreurs' ;
   'ERREUR' 5 ;
'FINSI' ;
'SI' interact ;
   'OPTION' 'DONN' 5 ;
'FINSI' ;   
*
* End of dgibi file couplage_TH1D_Th3D_2
*
'FIN' ;

 

 

© Cast3M 2003 - All rights reserved.
Disclaimer