Download Th1D-T3D-mono.dgibi

Back to the list

* fichier : Th1D-T3D-mono.dgibi
*
'OPTI' 'ECHO' 0 ;
**********************************************************************************************
*
* NOM DE FICHIER    : Th1D-T3D-mono.dgibi
*
* DESCRIPTION       : Modelisation thermique 3D - thermohydraulique 1D monophasique
*
**********************************************************************************************
* CAS TEST CORRESPONDANT AU RAPPORT DM2S/STMF/LMES/RT/12-013/A :
*                  Calculs thermohydrauliques HCLL/TBM realises avec le code CAST3M "fluide"
*                  Etude du refroidissement de la couverture tritigene du reacteur a fusion
*                  par un ecoulement d'helium
*
**********************************************************************************************
*       A U T E U R S
**********************************************************************************************
*
* Giacomo AIELLO    - CEA/DEN/DANS/DM2S/SEMT/BCCR - mel : giacomo.aiello@cea.fr
* Jean-Louis LASZLO - CEA/DEN/DANS/DM2S/STMF/LMES - mel : jean-louis.laszlo@cea.fr
* Stephane GOUNAND  - CEA/DEN/DANS/DM2S/SFMF/LMSF - mel : gounand@semt2.smts.cea.fr
*
**********************************************************************************************
**********************************************************************************************
*
*
* TEST RAPIDE       : la variable logique complet = FAUX;
* TEST NOMINAL      : la variable logique complet = VRAI;
*
*
**********************************************************************************************
*                   Exemple de resultats                                                     *
**********************************************************************************************
* Variable logique COMPLET                            FAUX              VRAI
*                  *******                            ****              ****
* Nombre iterations                                   1 seule           test
*                                                                       convergence
* Maillage                                            grossier          nominal
*
* Debit massique commun fluide 1 et 2                 0.11 kg/s         0.11 kg/s
* Flux impose                                         92.4 kW           92.4 kW
*
* Temperature minimale solide                        298.8 degres C    303.5 degres C
* Temperature maximale solide                        449.0 degres C    517.3 degres C
*
* Temperature minimale fluide 1                      300.0 degres C    300.0 degres C
* Temperature maximale fluide 1                      380.5 degres C    383.7 degres C
**********************************************************************************************
*
*
*                   VARIABLES LOGIQUES              
*
**********************************************************************************************

complet = FAUX ; comm TEST RAPIDE ;
*complet = VRAI ; comm TEST NOMINAL ; 

* EchColb = VRAI ;  comm loi echange thermique COLBURN ;
EchColb = FAUX ;  comm loi echange thermique GNIELINSKI ;


**********************************************************************************************
**********************************************************************************************
***                          PROCEDURES
**********************************************************************************************
**********************************************************************************************

debproc spabas spass*flottant ori*point ;
v9 = spass ** 0.5 ;
P2 = ori PLUS ( v9 0. 0. ) ;
Lig9 = DROI 1 ori P2 ;
sp9  = Lig9 TRAN 1 ( 0. v9 0. ) ;
elim sp9 1.E-6 ;
finpro sp9 ;

debproc loupe eye  c0*chpoint m0*maillage tjt0*mot gross0*flottant
        c/flottant ;
opti isov surf ;
* 1 pour UX - 2 pour UY - 3 pour UZ
xmt = 'COORDONNEE' 1  m0 ;
ymt = 'COORDONNEE' 2  m0 ;
dxmt = 'NOMC' 'UX' ('*' xmt gross0) 'NATURE' 'DISCRET' ;
dymt = 'NOMC' 'UY' ('*' ymt gross0) 'NATURE' 'DISCRET' ;
orig0 = 'FORME' ;
'FORME' dxmt ;
'FORME' dymt ;
SI ( EXIS c ) ;
 trac cach eye c0 m0 titr tjt0 ;
SINON ;
 trac cach eye c0 m0 (cont m0 ) titr tjt0 ;
FINSI ;
'FORME' orig0 ;
finproc ;

DEBP PARTSURF TABCONV*TABLE;
TABCONV . PART = TABLE;

NBOUCT = NBNO TABCONV.SURFCONV ;
REPE BOUC NBOUCT ;
    SI (&bouc EGA 1 );
        NEWPART = ELEM TABCONV . SURFCONV APPUYE
                            LARGEMENT TABCONV . SIN;
        TOTPART = NEWPART;
        PINI = BARY NEWPART;
    SINON ;
        NEWPART = ELEM TABCONV . SURFCONV APPUYE LARGEMENT TOTPART;
        SI ((NBNO NEWPART) EGA (NBNO TOTPART));
            QUITTER BOUC;
        SINON ;
            NEWPART = DIFF NEWPART TOTPART;
            TOTPART = TOTPART et NEWPART;
        FINSI;
    PFIN = BARY NEWPART;
    FINSI;
    TABCONV . PART . &bouc = NEWPART;
FIN bouc;
FINP;

DEBPROC FLUI1D_2 TABCONV*TABLE
      DECAX*FLOTTANT DECAY*FLOTTANT kdec*FLOTTANT spa9*FLOTTANT ;
TABCONV . FAC0 = TABLE ;
TABCONV . ELE0 = TABLE ;

DECAX = kdec * DECAX ;
DECAY = kdec * DECAY ;

Npart = DIME TABCONV .'PART' ;
PINI = BARY TABCONV .'PART'. 1 ;
PINI = PINI PLUS ( DECAX DECAY 0. ) ;
SurBas = spabas spa9 PINI ;
TABCONV .'BAFLU' = SurBas ;

*mess Npart ;

kfac = 0 ; comm initialisation compteur ;
REPE BOUC Npart ;
 NS = NBEL TABCONV .'PART'.&bouc ;
*mess NS ;
 Slat = ( mesu ( TABCONV .'PART'.&bouc ));
 HautF = Slat / ( 4. * ( spa9 ** 0.5 )) ;
* mess &bouc spa9 hautf slat ;
 opti elem cub8 ;
   si ( ega &bouc 1 ) ;
    NEWFLU = SurBas VOLU 1 'TRAN' ( 0. 0. HautF  ) ;
    TOTFLU = NEWFLU ;
    NEWSOL = TABCONV .'PART'.&bouc ;
    TOTSOL = NEWSOL ;
   SINON ;
    NEWFLU = ( FACE 2 NEWFLU ) volu 1 'TRAN' ( 0. 0. HautF  ) ;
    TOTFLU = TOTFLU ET NEWFLU ;
    NEWSOL = TABCONV .'PART'.&bouc ;
    TOTSOL = TOTSOL et NEWSOL ;
   FINSI ;

 repe w ( NBEL TABCONV .'PART'.&bouc ) ;
  kfac = kfac + 1 ;
* mess &bouc &w kfac ;
  TABCONV . FAC0 . kfac = TABCONV .'PART'.&bouc ELEM &w ;
  TABCONV . ELE0 . kfac = NEWFLU ;
 fin w ;

FIN bouc;

TABCONV .'MTFLU' = totflu ;
TABCONV .'OUFLU' = FACE 2 NEWFLU ;

FINP ;


debproc matcon0 tfa00*TABLE tel00*TABLE ;
*** CREATION DE LA MATRICE DE COUPLAGE                      *
*** ----------------------------------                      *
nitab = 'DIME' tfa00 ;
'REPETER' iitab nitab ;
   ifac = tfa00 . &iitab ;
   iqfac = 'CHANGER' ifac 'QUAF' ;
   miqfac = 'MODELISER' iqfac 'NAVIER_STOKES' discr ;
   ciqfac = ('DOMA' miqfac 'CENTRE') 'POIN' 1 ;
   iele = tel00 . &iitab ;
   iqele = 'CHANGER' iele 'QUAF' ;
   miqele = 'MODELISER' iqele 'NAVIER_STOKES' discr ;
   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) ;
      cfc00 = clfc ;
      mcf00 = mlcf ;
      mfc00 = mlfc ;
   'SINON' ;
      cfc00 = 'ET' cfc00 clfc ;
      mcf00 = 'ET' mcf00 mlcf ;
      mfc00 = 'ET' mfc00 mlfc ;
   'FINSI' ;
'FIN' iitab ;
finpro cfc00 mcf00 mfc00 ;

DEBPROC tcpu9 a/flottant ;  
tcpu = 'TEMPS' 'NOEC' ;    
tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ; 
SI ( exis a ) ;
 RESPRO tcpus ;
SINON ;
 'MESSAGE' ( et ('CHAINE' 'tcpus = ' tcpus) ' secondes') ;    
FINSI ;
FINPROC  ;  
DEBPROC PRPHE8 prt7C*'LISTREEL' roo*'LISTREEL' muo*'LISTREEL'
        cpo*'LISTREEL' lambo*'LISTREEL' T/'CHPOINT'
        TL/'LISTREEL'  TS/'FLOTTANT' ;
prt7 = prt7C + ( prog ( dime prt7C ) * 273.15 ) ; comm en Kelvin ;

   'SI' ('EXISTE' T) ;
     rhohe8 = 'IPOL' T prt7 roo ;
     muhe8  = 'IPOL' T prt7 muo ;
     lhe8   = 'IPOL' T prt7 lambo ;
     cphe8  = 'IPOL' T prt7 cpo ;
    'RESPRO' rhohe8 muhe8 lhe8 cphe8 ;
   'FINSI' ;

   'SI' ('EXISTE' TL) ;
     rhohe8 = 'IPOL' TL prt7 roo ;
     muhe8  = 'IPOL' TL prt7 muo ;
     lhe8   = 'IPOL' TL prt7 lambo ;
     cphe8  = 'IPOL' TL prt7 cpo ;
    'RESPRO' rhohe8 muhe8 lhe8 cphe8 ;
   'FINSI' ;

   'SI' ('EXISTE' TS) ;
     rhohe8 = 'IPOL' TS prt7 roo ;
     muhe8  = 'IPOL' TS prt7 muo ;
     lhe8   = 'IPOL' TS prt7 lambo ;
     cphe8  = 'IPOL' TS prt7 cpo ;
    'RESPRO' rhohe8 muhe8 lhe8 cphe8 ;
   'FINSI' ;
'FINPROC' ;

DEBPROC PTH_HE8 X/'CHPOINT' XL/'LISTREEL'  XS/'FLOTTANT' ;
TABHE8 = TABLE ;
   'SI' ('EXISTE' X) ;
    rhohe8 muhe8 lhe8 cphe8 = 
       PRPHE8 PROGTHe PROGRHO PROGMU PROGCP PROGLAM X ;
   'FINSI' ;

   'SI' ('EXISTE' XL) ;
    rhohe8 muhe8 lhe8 cphe8 = 
       PRPHE8 PROGTHe PROGRHO PROGMU PROGCP PROGLAM XL ;
   'FINSI' ;

   'SI' ('EXISTE' XS) ;
    rhohe8 muhe8 lhe8 cphe8 = 
       PRPHE8 PROGTHe PROGRHO PROGMU PROGCP PROGLAM XS ;
   'FINSI' ;
TABHE8 . 'ROG'  = rhohe8 ;
TABHE8 . 'MUG'  = muhe8 ;
TABHE8 . 'LAMG' = lhe8 ;
TABHE8 . 'CPG'  = cphe8 ;
FINPRO TABHE8 ;

DEBPROC PRACIER prt7C*'LISTREEL' roo*'LISTREEL' 
        cpo*'LISTREEL' lambo*'LISTREEL' T/'CHPOINT'
        TL/'LISTREEL'  TS/'FLOTTANT' ;
prt7 = prt7C + ( prog ( dime prt7C ) * 273.15 ) ; comm en Kelvin ;
   'SI' ('EXISTE' T) ;
     rhoaci = 'IPOL' T prt7 roo ;
     laci   = 'IPOL' T prt7 lambo ;
     cpaci  = 'IPOL' T prt7 cpo ;
    'RESPRO' rhoaci laci cpaci ;
   'FINSI' ;

   'SI' ('EXISTE' TL) ;
     rhoaci = 'IPOL' TL prt7 roo ;
     laci   = 'IPOL' TL prt7 lambo ;
     cpaci  = 'IPOL' TL prt7 cpo ;
    'RESPRO' rhoaci laci cpaci ;
   'FINSI' ;

   'SI' ('EXISTE' TS) ;
     rhoaci = 'IPOL' TS prt7 roo ;
     laci   = 'IPOL' TS prt7 lambo ;
     cpaci  = 'IPOL' TS prt7 cpo ;
    'RESPRO' rhoaci laci cpaci ;
   'FINSI' ;
'FINPROC' ;

DEBPROC PTH_ACI X/'CHPOINT' XL/'LISTREEL'  XS/'FLOTTANT' ;
TABACI = TABLE ;
  'SI' ('EXISTE' X) ;
    rhohe8 lhe8 cphe8 = 
       PRACIER PROGTAc ROACIER CPACIER LAMACIER X ;
   'FINSI' ;

   'SI' ('EXISTE' XL) ;
    rhohe8  lhe8 cphe8 = 
       PRACIER PROGTAc ROACIER CPACIER LAMACIER XL ;
   'FINSI' ;

   'SI' ('EXISTE' XS) ;
    rhohe8  lhe8 cphe8 = 
       PRACIER PROGTAc ROACIER CPACIER LAMACIER XS ;
   'FINSI' ;
TABACI . 'ROS'  = rhohe8 ;
TABACI . 'LAMS' = lhe8 ;
TABACI . 'CPS'  = cphe8 ;
FINPRO TABACI ;

'DEBPROC'  chpobor mess1/'MOT' T/'FLOTTANT' TC/'CHPOINT' ;
*Elle permet d'afficher les valeurs minimales et maximales d'un
*champoint avec message optionnel
MesEx0 = EXISTE mess1 ;
'SI' ('EXISTE' TC) ;
  SI MesEx0 ;
  mess mess1 ' : ' ( MINI TC ) ( MAXI TC) ;
 SINON ;
  mess ( MINI TC ) ( MAXI TC) ;
 FINSI ;
FINSi ;
'SI' ('EXISTE' T) ;
  SI MesEx0 ;
  mess mess1 ' : ' T ;
 SINON ;
  mess T ;
 FINSI ;
FINSi ;
FINPROC ;

debproc ecform c1*MOT c2*MOT c3*MOT v0*FLOTTANT ;
c1 = et c1 ' = ' ;
c3 = et ' ' c3 ;
c20 = ( chaine 'FORMAT' c2 v0) ;
mess (et (et c1 c20 ) c3) ;
finproc ;

**********************************************************************************************
**********************************************************************************************
***                       FIN  PROCEDURES
**********************************************************************************************
**********************************************************************************************


**********************************************************************************************
*
*                   DONNEES GEOMETRIE               
*
**********************************************************************************************

OPTION TRAC PSC ;

OPTI DIME 3;
OPTI ELEM TET4;

DENS 0.003;

R = 0.575; comm profondeur radiale ;
H = 1.660; comm hauteur totale ;

RFW = 0.100; comm rayon de courbure de la FW ;
PITCHFW = 0.0219; comm pitch entre canaux ;
EPFW = 0.030; comm epaisseur totale FW ;

DCT = 0.015; comm dimensions des canaux de refroidissement ;
DCP = 0.015;

DPL = 0.0035; comm distance des canaux de la surface exposee au plasma ;

DCH = (PITCHFW - DCT) / 2.;


EPS0 = 1.E-6 ;

***********************************************************************
*
*                   MAILLAGE               
*
***********************************************************************

si complet ;
 N1V1 = 3 ;
 N2V1 = 20 ;
 NV2  = 8 ;
 NV3  = 40 ;
 NV7  = 6 ;
sinon ;
* N1V1 = 3 ;
* N2V1 = 3 ;
* NV2  = 4 ;
* NV3  = 4 ;
* NV7  = 6 ;
* modif gounand
 N1V1 = 2 ;
 N2V1 = 2 ;
 NV2  = 2 ;
 NV3  = 2 ;
 NV7  = 2 ;
FINSI ;

NV4  = NV2  ;
NV5  = N2V1 ;
NV6  =  N1V1 ;
NV8  =  N1V1 ;
NV9  = N2V1 ;
NV10 = NV2 ;
NV11 = NV3 ;
NV12 = NV2 ;
N1V13 = N2V1 ;
N2V13 = N1V1 ;


P1 = 0. 0. 0. ;
P2 = 0. PITCHFW 0.;
P3 = 0. PITCHFW ( -1. * EPFW );
P4 = 0. 0. (-1. * EPFW);


PC1 = 0. DCH ( -1. * DPL);
PC2 = 0. (DCH + DCT) ( -1. * DPL);
PC3 = 0. (DCH + DCT) ( -1. * (DPL + DCP) );
PC4 = 0. DCH ( -1. * (DPL + DCP) );




L1 = D P1 P2;
L2 = D P2 P3;
L3 = D P3 P4;
L4 = D P4 P1;

LC1 = D PC1 PC2;
LC2 = D PC2 PC3;
LC3 = D PC3 PC4;
LC4 = D PC4 PC1;

C1 = L1 et L2 et L3 et L4;
C2 = LC1 et LC2 et LC3 et LC4;

S1 = SURF ( C1 et (INVE C2) ) PLAN;
V1 = VOLU TRAN N1V1 (PITCHFW 0. 0.) S1;
V1 = VOLU TRAN N2V1 ( (R - RFW - PITCHFW) 0. 0.) V1;


S2 = FACE 2 V1;
V2 = VOLU ROTA NV2 90. S2 
        ( (R - RFW) 0. (-1 * RFW))
        ( (R - RFW) PITCHFW (-1 * RFW));



S3 = FACE 2 V2;
V3 = VOLU TRAN NV3 ( 0. 0. (-1 * (H - (2. * RFW)))) S3;

S4 = FACE 2 V3;
V4 = VOLU ROTA NV4 90. S4 
        ( (R - RFW) 0. (-1 * (H - RFW)))
        ( (R - RFW) PITCHFW (-1 * (H - RFW)));

S5 = FACE 2 V4;
V5 = VOLU TRAN NV5 ( (-1. * (R - RFW - PITCHFW) ) 0. 0.) S5;


S6 = FACE 2 V5;
S7 = PROJ (1. 0. 0.) S6 PLAN 
                ( (0. 0. (-1. * H) ))  
                ( (PITCHFW) (PITCHFW) (-1. * H) )  
                ( 0. 0. 0.);  
V6 = VOLU NV6 S6 S7;


S8 = PROJ S7 (0. 1. 0.) PLAN
                ( 0. (2*PITCHFW) H )  
                ( PITCHFW (PITCHFW) H )  
                ( 0. (2*PITCHFW) 0.);  
V7 = VOLU NV7 S7 S8;



S9 = PROJ S8 (-1. 0. 0.) PLAN
                    (PITCHFW. 0. 0.)
                    (PITCHFW 1. 0.)
                    (PITCHFW 0. 1.);
V8 = VOLU NV8 S8 S9;
                    
V9 = VOLU TRAN NV9 ( (R - RFW - PITCHFW)  0. 0.) S9; 

S10 = FACE 2 V9;
V10 = VOLU ROTA NV10 -90. S10 
        ( (R - RFW) 0. (-1 * (H - RFW)))
        ( (R - RFW) PITCHFW (-1 * (H - RFW)));

S11 = FACE 2 V10;
V11 = VOLU TRAN NV11 ( 0. 0.  (H - (2. * RFW))) S11;


S12 = FACE 2 V11;
V12 = VOLU ROTA NV12 -90. S12 
        ( (R - RFW) 0. (-1 * RFW))
        ( (R - RFW) PITCHFW (-1 * RFW));

S13 = FACE 2 V12;
V13 = VOLU TRAN N1V13 ( (-1. * (R - RFW - PITCHFW)) 0. 0.) S13;
V13 = VOLU TRAN N2V13 ( (-1.*PITCHFW) 0. 0.) V13;


VTUBE = V1 ET V2 ET V3 ET V4 ET V5 ET V6 ET V7 ET V8
        ET V9 ET V10 ET V11 ET V12 ET V13;
        
ELIM EPS0 VTUBE;

mess ( nbel vtube ) ( nbno vtube );


VTUBE = REGE VTUBE;

VTUBE2 = VTUBE SYME PLAN 
        (0. 0. (H / -2.))
        (0. 1. (H / -2.))
        (1. 0. (H / -2.));
            
VTUBE2 = DEPL PLUS VTUBE2 (0. (2*PITCHFW) 0.);

VFW = VTUBE et VTUBE2;
ELIM EPS0 VFW;
VFW = REGE VFW;

TRAC CACH VFW;


OPTI SAUV 'MAILLAGE.SAUV';
SAUV VFW;

*******************************************************************
*******************************************************************
***                                                              **
***      COUPLAGE THERMIQUE 3D - THERMOHYDRAULIQUE 1D            **
***                                                              **
*******************************************************************
*******************************************************************

OPTI DIME 3 ELEM CUB8 ;


OPTI REST 'MAILLAGE.SAUV'; comm Restitution du maillage ;
REST;


*** #################################################################*
*** ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*
***                          DECLARATIONS
*** ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*
*** #################################################################*


oeilx = 1.e30 0.    0.    ;
oeily = 0.    1.e30 0.    ;
oeilz = 0.    0.    1.e30 ;
oeil = 1.e30 1.e30 1.e30 ;
epselim = 1.e-7 ;

discr   = 'LINE' ; comm elements lineaires ;
* discr = 'QUAF' ; comm elements quadratiques ;

nomsov='Th1D-T3D-mono.sauv' ;

TFLUec  = 300.   ;
PsatB  = 1.75  ; comm pression absolue en bars ;
ppa   = PsatB * 1.e5 ; comm pression absolue en Pa ;
CK     = 273.15 ; comm Celsius en Kelvin ;
TFLUek   = TFLUec + CK ;
F_SPLSM = 0.5E6 ; comm flux surface plasma ;
F_SBZ   = 0.1E6 ; comm flux surface interne ;

nitr1 = 1 ;
omgr1 = 1. ;
SI (non complet)  ;
 nb_bouc = 1 ;
SINON ;
 nb_bouc = 40 ; comm boucle de calcul ;
FINSI ;

Nusslam = 4.364 ; comm Nusselt en laminaire tube circulaire ;
seuilam = 2300. ; comm transition laminaire turbulent  ;

*** Loi de convection forcee en regime turbulent (COLBURN)
A_colb = 0.023 ;
B_colb = 0.8D0 ;
C_colb = ('/' 1.D0 3.D0) ;

*** Loi de convection forcee en regime turbulent (Gnielinski)
A_gnie = 1.7372 ;
B_gnie = 1.964 ;
C_gnie = -3.8215 ;
s_23 = 2. / 3. ;


*** Methode d'inversion
MetInv = 'DIRECT' ; comm robuste mais lent ;
*** MetInv = 'ITERAT' ; comm moins robuste mais plus rapide ;
SI ( EGA MetInv 'DIRECT' ) ;
 IndInv = 1 ;
SINON ;
 IndInv = 3 ;
FINSI ;

RegEcoM = TABLE ;
RegEcoD = TABLE ;
RTUR = 'TURB' ;
RLAM = 'LAMI' ;


***     DONNEES RELATIVES A L ECOULEMENT:
GTUBE = 1.3 / 12. ; comm debit massique kg par seconde ;
SEC = 0.015 * 0.015 ; comm section de passage ;
PMOUL = 0.015 * 4. ; comm perimetre mouille ;
DIAM = 4 * SEC / PMOUL ; comm diametre hydraulique ;

*****************************************************************
***        PROPRIETES MATERIAUX
*****************************************************************
***       ici les temperatures PROGTAc sont en degres CELSIUS
PROGTAc  = PROG 20.  100. 200. 300. 400. 500. 600. ;
LAMACIER = PROG 25.9 27.0 28.1 28.8 29.2 29.0 28.5 ;
CPACIER  = PROG 446.  486.  524.  564.  621. 701.  805. ;
ROACIER  = PROG 7758. 7756. 7752. 7750. 7746.7 7743. 7740. ;


*** Conducibilite Thermique ACIER
EVOCDST    = EVOL MANU
 'T' PROGTAc
 'K' LAMACIER ;

*** On n'utilise pas ces valeurs
*** Chaleur specifique ACIER
EVOCPST   = EVOL MANU
 'T' PROGTAc
 'C' CPACIER  ;
*** DENSITE ACIER
EVOROST    = EVOL MANU
'T'   PROGTAc
'RHO' ROACIER  ;


*
***                            HELIUM 8 MPA
*
*** ATTENTION: LES PROGRESSIONS DOIVENT AVOIR LA MEME LONGEUR
***       ici les temperatures PROGTHe sont en degres CELSIUS
PROGTHe = PROG 250. 260. 270. 280. 290. 300. 310. 320. 330. 340.
              350. 360. 370. 380. 390. 400. 410. 420. 430. 440. 450.
              460. 470. 480. 490. 500. 510. 520. 530. 540. 550. ;
PROGRHO = PROG
7.2195 7.0872 6.9596 6.8365 6.7176 6.6029 6.4919 6.3847 6.2809 6.1804
6.0831 5.9888 5.8974 5.8087 5.7227 5.6391 5.5580 5.4791 5.4025 5.3280
5.2554 5.1849 5.1162 5.0493 4.9841 4.9206 4.8587 4.7983 4.7394 4.6819
4.6258 ;

PROGCP = PROG
 5183.9 5184.2 5184.5 5184.8 5185.1 5185.4 5185.6 5185.9 5186.1 5186.3
 5186.5 5186.7 5186.9 5187.1 5187.3 5187.4 5187.6 5187.7 5187.9 5188.0
 5188.1 5188.3 5188.4 5188.5 5188.6 5188.7 5188.8 5188.9 5188.9 5189.0
 5189.1  ;

PROGLAM = PROG
  0.23121 0.23432 0.23741 0.24049 0.24354
  0.24658 0.24960 0.25261 0.25560 0.25857
  0.26152 0.26446 0.26738 0.27029 0.27318
  0.27606 0.27892 0.28176 0.28460 0.28741
  0.29022 0.29300 0.29578 0.29854 0.30129
  0.30402 0.30675 0.30946 0.31215 0.31484 0.31751;

PROGMU = PROG
  2.86840E-05  2.90485E-05  2.94108E-05  2.97711E-05  3.01293E-05
  3.04856E-05  3.08401E-05  3.11927E-05  3.15435E-05  3.18927E-05
  3.22402E-05  3.25862E-05  3.29307E-05  3.32737E-05  3.36154E-05
  3.39557E-05  3.42947E-05  3.46325E-05  3.49691E-05  3.53045E-05
  3.56389E-05  3.59722E-05  3.63046E-05  3.66360E-05  3.69664E-05
  3.72960E-05  3.76248E-05  3.79527E-05  3.82799E-05  3.86064E-05
  3.89322E-05;

EVOROHE  = EVOL MANU 'T' PROGTHe 'RHO' PROGRHO;
EVOCPHE  = EVOL MANU 'T' PROGTHe 'CP' PROGCP;
EVOMUHE  = EVOL MANU 'T' PROGTHe 'MU' PROGMU;
EVOLAMHE = EVOL MANU 'T' PROGTHe 'LAM' PROGLAM;


*** helium : la procedure PTH_HE8 sort ro mu lambda cp
*** acier  : la procedure PTH_ACI sort ro lambda cp


****************************************************************
***       IDENTIFICATION DES SURFACES D ECHANGE
****************************************************************


XMIN = MINI (COOR 1 VFW);
ELXMIN = ELEM (ENVE VFW) APPUYE
                (POINT VFW PLAN
                        (XMIN 0. 0.)
                        (XMIN 1. 0.)
                        (XMIN 0. 1.));

CTUBES = CCON (CONT ELXMIN);
EXPTUBES = CCON (DIFF (ENVE VFW) ELXMIN);

list EXPTUBES ;


SIN1 = CTUBES . 2;
SIN2 = CTUBES . 6;
SCONV1 = EXPTUBES . 2;
SCONV2 = EXPTUBES . 3;

****************************************************************
***       PARTITIONNEMENT DES SURFACES D ECHANGE
****************************************************************


TUBE1 = TABLE;
TUBE2 = TABLE;

TUBE1 . 'G' = GTUBE;
TUBE1 . 'PM' = PMOUL;
TUBE1 . 'DH' = DIAM;
TUBE1 . 'RHO' = EVOROHE;
TUBE1 . 'LAM' = EVOLAMHE;
TUBE1 . 'CP' = EVOCPHE;
TUBE1 . 'MU' = EVOMUHE;

TUBE2 . 'G' = GTUBE;
TUBE2 . 'PM' = PMOUL;
TUBE2 . 'DH' = DIAM;
TUBE2 . 'RHO' = EVOROHE;
TUBE2 . 'LAM' = EVOLAMHE;
TUBE2 . 'CP' = EVOCPHE;
TUBE2 . 'MU' = EVOMUHE;

Spass1 = TUBE1 . 'PM' * TUBE1 . 'DH' / 4. ;
Spass2 = TUBE2 . 'PM' * TUBE2 . 'DH' / 4. ;


TUBE1 . SIN = SIN1;
TUBE2 . SIN = SIN2;
TUBE1 . SURFCONV = SCONV1;
TUBE2 . SURFCONV = SCONV2;

*** Flux de chaleur sur la FW
YMAX = MAXI (COOR 2 VFW);
YMIN = MINI (COOR 2 VFW);

ELYMIN = ELEM (ENVE VFW) APPUYE
            (POINT VFW PLAN
                (0. YMIN 0.)
                (1. YMIN 0.)
                (0. YMIN 1.) );

ELYMAX = ELEM (ENVE VFW) APPUYE
            (POINT VFW PLAN
                (0. YMAX 0.)
                (1. YMAX 0.)
                (0. YMAX 1.) );

ENVEXT = DIFF (ENVE VFW)
             (ELYMAX et ELYMIN et ELXMIN et
              TUBE1 . SURFCONV et TUBE2 . SURFCONV);

SFW2 = (CCON ENVEXT) . 2;
SFW1 = (CCON ENVEXT) . 1;

ZMAX = MAXI (COOR 3 SFW2);
ELSIDE = ELEM SFW2 APPUYE
            (POINT SFW2 PLAN
                (0. 0. ZMAX)
                (1. 0. ZMAX)
                (0. 1. ZMAX) 0.0001);

XMAX = MAXI (COOR 1 ELSIDE);
SPLSM = ELEM SFW2 APPUYE
            (POINT (COOR 1 SFW2) EGSUPE (XMAX - 0.0001));

xmil ymil zmil = coor splsm ;
pmil = ( maxi xmil ) (0.5 * ((mini ymil) + ( maxi ymil)))
  (0.5 * ((mini zmil ) + ( maxi zmil))) ;

ZMAX = MAXI (COOR 3 SFW1);
ELSIDE = ELEM SFW1 APPUYE
            (POINT SFW1 PLAN
                (0. 0. ZMAX)
                (1. 0. ZMAX)
                (0. 1. ZMAX) 0.0001);

XMAX = MAXI (COOR 1 ELSIDE);
SBZ = ELEM SFW1 APPUYE
            (POINT (COOR 1 SFW1) EGSUPE (XMAX - 0.0001));

W_SPLSM = F_SPLSM * ( mesu SPLSM ) ;
W_SBZ   = F_SBZ   * ( mesu SBZ ) ;

XTT YTT ZTT = COOR VFW ;
decx = 2. * (( MAXI XTT ) - ( MINI XTT )) ;
decy = 2. * (( MAXI YTT ) - ( MINI YTT )) ;

PARTSURF TUBE1 ;
FLUI1D_2 TUBE1 decx decy 1.0 Spass1 ;

PARTSURF TUBE2;
FLUI1D_2 TUBE2 decx decy 1.1 Spass2 ;



sc1 = tube1.surfconv ;
sc2 = tube2.surfconv ;

*trac cach ( vfw et tube1.mtflu et tube2.mtflu ) ;


elim epselim  ( tube1.baflu et tube1.mtflu et sc1 ) ;
pelim1 = '*' ('**' ('/' ('MESURE' tube1.mtflu) ('NBEL' tube1.mtflu))
                  ('/' 1.D0 ('VALEUR' 'DIME'))) 1.D-3 ;

'ELIMINATION' (sc1 et VFW ) pelim1 ;
'ELIMINATION' (tube1.baflu et tube1.mtflu  ) pelim1 ;

elim epselim  ( tube2.baflu et tube2.mtflu et sc2 ) ;
pelim2 = '*' ('**' ('/' ('MESURE' tube2.mtflu) ('NBEL' tube2.mtflu))
                  ('/' 1.D0 ('VALEUR' 'DIME'))) 1.D-3 ;

'ELIMINATION' (sc2 et VFW ) pelim2 ;
'ELIMINATION' (tube2.baflu et tube2.mtflu  ) pelim2 ;

*** On va construire le chpoint des produits scalaires de la normale
*** a la face SPLSM et le vecteur unitaire i (1. 0. 0. ), egal au cosinus
*** de la normale et du vecteur i, pour ponderer le flux impose et 
*** faire en sorte que sa direction soit donc constante. On cree un
*** nouveau maillage SPLSM2 oriente vers l'exterieur. _SPLSM2 et $SPLSM2
*** sont respectivement les elements QUAF et le modele Navier-Stokes
*** correspondant.
*** Idem avec la face SBZ

*** On va multiplier ce produit scalaire par le flux impose.

ListXm = prog ;
repe x ( nbno splsm ) ;
 ListXm = ListXm et ( prog 1. 0. 0. ) ;
fin x ;
splsm = inve ( orie  SPLSM 'POINT' ( 0. 0. 0. )) ;
SBZ   = inve ( orie  SBZ   'POINT' ( 0. 0. 0. )) ;
CHPx1 = MANU 'CHPO' splsm 
          ('MOTS' 'UX' 'UY' 'UZ') ListXm ;
CHPx2 = MANU 'CHPO' sbz 
          ('MOTS' 'UX' 'UY' 'UZ') ListXm ;

*** Generation des QUAFs
*** --------------------
_baflu1     = 'CHANGER' tube1.baflu   'QUAF' ;
_mtflu1     = 'CHANGER' tube1.mtflu   'QUAF' ;
_Clim1      = 'CHANGER' sc1           'QUAF' ;
_VFW        = 'CHANGER' VFW           'QUAF' ;

_baflu2     = 'CHANGER' tube2.baflu   'QUAF' ;
_mtflu2     = 'CHANGER' tube2.mtflu   'QUAF' ;
_Clim2      = 'CHANGER' sc2           'QUAF' ;
_SPLSM      = 'CHANGER' SPLSM         'QUAF' ;
_SBZ        = 'CHANGER' SBZ           'QUAF' ;


'TEMPS' 'ZERO' ;
cfc1 mcf1 mfc1 = matcon0 ( tube1 . fac0 ) ( tube1 . ele0 ) ;

*** Elimination des points doubles
*** ------------------------------
'ELIMINATION' ( _baflu1  'ET' _mtflu1  'ET' cfc1 ) pelim1 ;
'ELIMINATION' ( _Clim1 'ET' _VFW 'ET' cfc1 ) pelim1 ;

mcf1 = 'CHANGER' mcf1 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc1 = 'CHANGER' mfc1 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;

cfc2 mcf2 mfc2 = matcon0 ( tube2 . fac0 ) ( tube2 . ele0 ) ;




*** Elimination des points doubles
*** ------------------------------
'ELIMINATION' ( _baflu2  'ET' _mtflu2  'ET' cfc2 ) pelim2 ;
'ELIMINATION' ( _Clim2 'ET' _VFW 'ET' cfc2 ) pelim2 ;

mcf2 = 'CHANGER' mcf2 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc2 = 'CHANGER' mfc2 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
TCpu1 = tcpu9 0. ; comm temps ecoule ;



*** ***** Densite de puissance neutronique
afw = 2.7116 ;
bfw = 13.673 ;
cfw = 3.3047 ;
dfw = 5.425  ;
FACTF = 0.78/0.75;


Pr_Xfw = PROG  0. PAS 0.001 1 ;
Pr_NWLfw  = FACTF * (  (afw * (EXP ( -1. * bfw * Pr_Xfw))) +
                       (cfw * (EXP ( -1. * dfw * Pr_Xfw)))
                    );
EVQfw = EVOL MANU 'XRAD' Pr_Xfw 'FW' Pr_NWLfw  ;
X_TOT = (COOR 1 VFW);
NWP_fw = IPOL X_TOT Pr_Xfw Pr_NWLfw;
NWP_fw = NWP_fw * 1.E6;



*** Recapitulatif
*** -------------
opti echo 0 ;
Mess 'Recapitulatif : ' ;
Mess '-------------' ;
Mess ;

SI ( EGA MetInv 'DIRECT' ) ;
 Mess 'Methode d inversion directe, robuste mais lente' ;
SINON ;
 Mess 
'Methode d inversion iterative, moins robuste mais plus rapide';
FINSI ;

SI EchColb ;
 MESS 'Loi echange turbulent COLBURN' ;
SINON ;
 MESS 'Loi echange turbulent Gnielinski' ;
FINSI ;

SI (non complet) ;
 MESS 'Calcul test avec ' nb_bouc ' iterations' ;
SINON ;
 MESS 'Calcul nominal avec ' nb_bouc ' iterations' ;
FINSI ;

*gounand opti echo 1 ;

***********************************************************
*                                                         *
*** PREPARATION DU CALCUL                                 *
*                                                         *
***********************************************************
*
*** Modeles Navier-Stokes
*** ---------------------

$baflu1   =  'MODELISER' _baflu1  'NAVIER_STOKES' discr  ;
$mtflu1   =  'MODELISER' _mtflu1  'NAVIER_STOKES' discr  ;
$Clim1    =  'MODELISER' _Clim1   'NAVIER_STOKES' discr  ;
$baflu2   =  'MODELISER' _baflu2  'NAVIER_STOKES' discr  ;
$mtflu2   =  'MODELISER' _mtflu2  'NAVIER_STOKES' discr  ;
$Clim2    =  'MODELISER' _Clim2   'NAVIER_STOKES' discr  ;
$VFW      =  'MODELISER' _VFW     'NAVIER_STOKES' discr  ;
$SPLSM    =  'MODELISER' _SPLSM   'NAVIER_STOKES' discr  ;
$SBZ      =  'MODELISER' _SBZ     'NAVIER_STOKES' discr  ;


FSPLSM9 = F_SPLSM * ( PSCA  ( doma $SPLSM NORMALEV ) 
CHPx1 ('MOTS' 'UX' 'UY' 'UZ') ('MOTS' 'UX' 'UY' 'UZ')) ;
FSBZ9 = F_SBZ * ( PSCA  ( doma $SBZ NORMALEV ) 
CHPx2 ('MOTS' 'UX' 'UY' 'UZ') ('MOTS' 'UX' 'UY' 'UZ')) ;

trac cach FSPLSM9 splsm titr 'flux impose sur face SPLSM - W/m2' ;
trac cach FSBZ9   sbz   titr 'flux impose sur face  SBZ - W/m2' ;


NWP_fwe = NOEL $VFW NWP_fw ; comm sources ;
VoSol = ( maxi ( resu ( DOMA $VFW 'VOLUME' ))) ;
SouSol = ( maxi ( resu  NWP_fw  )) ;



*** Initialisation des tables pour les calculs Thermiques structures et fluides

rvsol = 'EQEX' 'ITMA' 1 'NITER' nitr1 'FIDT' 1000 'OMEGA' omgr1 ;

rvsol = 'EQEX' rvsol
         'OPTI' 'EF' 'IMPL'
      'ZONE' $VFW     'OPER' 'LAPN' 'Lacier'   'INCO' 'TSOL'
      'ZONE' $Clim1   'OPER' 'ECHI' 'H1' 'T1'  'INCO' 'TSOL' 'TSOL'
      'ZONE' $Clim2   'OPER' 'ECHI' 'H2' 'T2'  'INCO' 'TSOL' 'TSOL'
        ;


rvsol = 'EQEX' rvsol
         'OPTI' 'EF' 'IMPL'
         'ZONE' $VFW  'OPER' 'FIMP' 'Soursol' 'INCO' 'TSOL'
         'ZONE' $SPLSM 'OPER' 'FIMP' FSPLSM9  'INCO' 'TSOL'
         'ZONE' $SBZ   'OPER' 'FIMP' FSBZ9    'INCO' 'TSOL'  
        ;

rvsol . 'METHINV' . 'TYPINV' = IndInv ;

*** Dans l'eau

rvflu = 'EQEX' 'ITMA' 1 'NITER' 1 'FIDT' 1000  ;
rvflu = 'EQEX' rvflu
         'OPTI' 'EF' 'IMPL' 'CENTREE' 'NOCONS'
'ZONE' $mtflu1 'OPER' 'KONV' 'RCPFLU1' 'UFLU1' 'LFLU1' 'INCO' 'TFLU'
'ZONE' $mtflu2 'OPER' 'KONV' 'RCPFLU2' 'UFLU2' 'LFLU2' 'INCO' 'TFLU'
        ;

rvflu = 'EQEX' rvflu
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtflu1 'OPER' 'LAPN' 'LFLU1' 'INCO' 'TFLU'
         'ZONE' $mtflu1 'OPER' 'FIMP' 'PFLU1' 'INCO' 'TFLU'
         'ZONE' $mtflu2 'OPER' 'LAPN' 'LFLU2' 'INCO' 'TFLU'
         'ZONE' $mtflu2 'OPER' 'FIMP' 'PFLU2' 'INCO' 'TFLU'
         ;


*** Conditions aux limites en temperature lames d'eau

rvflu = 'EQEX' rvflu
         'CLIM' 'TFLU' 'TIMP' ('DOMA' $baflu1 'MAILLAGE') TFLUek
         'CLIM' 'TFLU' 'TIMP' ('DOMA' $baflu2 'MAILLAGE') TFLUek
         ;

*** Initialisation de toutes les variables avant le calcul

rvi = 'TABLE' 'INCO' ;
rvsol . 'INCO' = rvi ;
rvflu . 'INCO' = rvi ;
rvi . 'TSOL' = 'KCHT' $VFW 'SCAL' 'SOMMET' TFLUek ;

LamAci = ( pth_aci TFLUek ) . 'LAMS' ;

rvi . 'Lacier'    = 'KCHT' $VFW   'SCAL' 'CENTRE' LamAci  ;
rvi . 'Soursol'  = 'KCHT' $VFW   'SCAL' 'CENTRE' NWP_fwe ;


rvi . 'T1'  = 'KCHT'  $Clim1 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F1'  = 'KCHT'  $Clim1 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H1'  = 'KCHT'  $Clim1 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T2'  = 'KCHT'  $Clim2 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F2'  = 'KCHT'  $Clim2 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H2'  = 'KCHT'  $Clim2 'SCAL' 'CENTRE' 0.D0 ;

rvi . 'TFLU' = 'KCHT' ( $mtflu1 et $mtflu2 ) 'SCAL' 'SOMMET' TFLUek ;
rvi . 'PFLU1'= 'KCHT' $mtflu1 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PFLU2'= 'KCHT' $mtflu2 'SCAL' 'CENTRE' 0.D0 ;

*** Copy of the fields to do the convergence check :
thep = 'COPIER' (rvi . 'TFLU') ;
tcbp = 'COPIER' (rvi . 'TSOL') ;

*** Initialisation des proprietes physiques de l'helium :
Rt_He8 = PTH_HE8 TFLUek ;
rflu  = Rt_He8 . 'ROG' ;
muflu = Rt_He8 . 'MUG' ;
lflu  = Rt_He8 . 'LAMG' ;
cpflu = Rt_He8 . 'CPG' ;

rvi . 'LFLU1'   =  'KCHT' $mtflu1 'SCAL' 'CENTRE' lflu ;
rvi . 'RCPFLU1' =  'KCHT' $mtflu1 'SCAL' 'CENTRE' ('*' rflu cpflu) ;
rvi . 'LFLU2'   =  'KCHT' $mtflu2 'SCAL' 'CENTRE' lflu ;
rvi . 'RCPFLU2' =  'KCHT' $mtflu2 'SCAL' 'CENTRE' ('*' rflu cpflu) ;


*** TUBEx . 'G' est un debit massique en kg/s
uflu1 = (( TUBE1 . 'G'  / Spass1 ) / rflu) ;
uflu2 = (( TUBE2 . 'G'  / Spass2 ) / rflu) ;
rvi . 'UFLU1' = 'KCHT' $mtflu1 'VECT' 'SOMMET' ( 0.D0 0.D0 uflu1 ) ;
rvi . 'UFLU2' = 'KCHT' $mtflu2 'VECT' 'SOMMET' ( 0.D0 0.D0 uflu2 ) ;

reyn6 = rflu * uflu1 * diam / muflu ;
mess diam reyn6  uflu1 ;




****************************************************************************
*** BOUCLE DE CALCUL                                                       *
***                                                                        *
****************************************************************************

'TEMPS' 'ZERO' ;

'REPETER' bouc nb_bouc ;

MESS' 1) Calcul de la Tre moy dans les fluides' ;
   tmyflu1 = 'NOEL' $mtflu1 (rvi . 'TFLU') ;
   t_flu1 = ( REDU (rvi . 'TFLU') tube1.mtflu ) ;
   tmyflu2 = 'NOEL' $mtflu2 (rvi . 'TFLU') ;
   t_flu2 = ( REDU (rvi . 'TFLU') tube2.mtflu ) ;

   'MESSAGE' 'Maxi Tre moy fluide 1 (degres C)     : '
             ('-' ('MAXIMUM' tmyflu1) CK) ;
   'MESSAGE' 'Maxi Tre moy fluide 2 (degres C)     : '
             ('-' ('MAXIMUM' tmyflu2) CK) ;
   chpobor 'Tre fluide 1 mini & maxi en C    ' ( t_flu1 - CK );
   chpobor 'Tre fluide 2 mini & maxi en C    ' ( t_flu2 - CK );


MESS '*  2) Transfert sur l interface';
   tmtif1  = '*' mcf1 tmyflu1 ;
   tmtif2  = '*' mcf2 tmyflu2 ;
   rvi . 'T1'  = 'KCHT' $clim1 'SCAL' 'CENTRE' tmtif1 ;
   rvi . 'T2'  = 'KCHT' $clim2 'SCAL' 'CENTRE' tmtif2 ;
chpobor 'Tre interface 1 (Celsius )' (( rvi . 'T1' ) - CK ) ;
chpobor 'Tre interface 2 (Celsius )' (( rvi . 'T2' ) - CK ) ;

MESS ' 3) Calcul du coefficient d echange a l interface   ' ;
***   3) Calcul du coefficient d'echange a l'interface
***      et des proprietes physiques du solide
***      a partir des temperatures a l'instant precedent.
*** Temperature de melange
   tm1 = rvi . 'T1' ;
   tm2 = rvi . 'T2' ;
*** Temperature a la paroi
   tp1 = 'NOEL' $clim1   (rvi . 'TSOL') ;
   tp2 = 'NOEL' $clim2   (rvi . 'TSOL') ;
*** Temperature de film
   tfi1 = '*' ('+' tm1 tp1) 0.5D0 ;
   tfi2 = '*' ('+' tm2 tp2) 0.5D0 ;   
   pptfil1 = PTH_HE8 tfi1 ;
   mufil1 = pptfil1 . 'MUG' ;
   lfil1  = pptfil1 . 'LAMG';
   cpfil1 = pptfil1 . 'CPG' ;
   pptfil2 = PTH_HE8 tfi2 ;
   mufil2 = pptfil2 . 'MUG' ;
   lfil2  = pptfil2 . 'LAMG';
   cpfil2 = pptfil2 . 'CPG' ;
   pranfil1 = '*' ('*' mufil1 cpfil1) ('INVERSE' lfil1) ;
   pranfil2 = '*' ('*' mufil2 cpfil2) ('INVERSE' lfil2) ;


SI ( EGA &bouc 1 ) ;
*** Calcul des nombres de Reynolds a partir de la temperature fluide (une fois)
opti echo  0 ;

   pptfil9 = PTH_HE8 ( mini t_flu1 ) ;
   rofil9 = pptfil9 . 'ROG' ; 
   mufil9 = pptfil9 . 'MUG' ;
   lfil9  = pptfil9 . 'LAMG';
   cpfil9 = pptfil9 . 'CPG' ;
   Reflu1 = ( ( TUBE1 . 'G' ) * ( TUBE1 . 'DH' ) / Spass1 ) / mufil9 ;

   pptfil9 = PTH_HE8 ( mini t_flu2 ) ;
   rofil9 = pptfil9 . 'ROG' ; 
   mufil9 = pptfil9 . 'MUG' ;
   lfil9  = pptfil9 . 'LAMG';
   cpfil9 = pptfil9 . 'CPG' ;
   Reflu2 = ( ( TUBE2 . 'G' ) * ( TUBE2 . 'DH' ) / Spass2 ) / mufil9 ;
Mess 'Nb de Reynolds ' ;
chpobor 'interne fluide 1 ' Reflu1 ; 
chpobor 'interne fluide 2 ' Reflu2 ; 

SI ( Reflu1 > seuilam ) ;
 RegEcoM.1 = RTUR ;
SINON ;
 RegEcoM.1 = RLAM ;
FINSI ;

SI ( Reflu2 > seuilam ) ;
 RegEcoM.2 = RTUR ;
SINON ;
 RegEcoM.2 = RLAM ;
FINSI ;

MESS ;
*gounand opti echo 1 ;
FINSI ;  comm Calcul des nombres de Reynolds ;

SI ( EGA RegEcoM.1 RTUR ) ; comm fluide UN  ;

SI EchColb ;
*** Colburn
 nuss1 = A_colb '*' ('**' Reflu1 B_colb)    
             '*' ('**' pranfil1 C_colb) ;    
SINON ;

*** Gnielinski
 ReyGnie  = Reflu1 - 1000. ;
 ReyLogG  = LOG Reflu1 ;

 K_gnie = A_gnie * (LOG (Reflu1 / 
               (( B_gnie * ReyLogG ) + C_gnie ))) ;
 F_frot = 0.5 / ( K_gnie ** 2. ) ;
 F_frot2 = F_frot ** 0.5 ;

 chpobor 'Reyn - 1000. ' ReyGnie ;
 chpobor 'coeff frott  ' ( 2. * F_frot );

 PrGnie = ((pranfil1 ** s_23 ) - 1. ) ;
 chpobor 'Pr**2/3 - 1  ' PrGnie ;
 nuss1 = F_frot * pranfil1 * ReyGnie / 
            ( 1. + (12.7 * F_frot2 * PrGnie )) ;
FINSI ;

SINON ;
 nuss1 = Nusslam ;
 nuss2 = Nusslam ;
FINSI ;

SI ( EGA RegEcoM.1 RTUR ) ; comm fluide DEUX  ;

SI EchColb ;
*** Colburn
 nuss2 = A_colb '*' ('**' Reflu2 B_colb)    
             '*' ('**' pranfil2 C_colb) ;   
SINON ;

*** Gnielinski
 ReyGnie  = Reflu2 - 1000. ;
 ReyLogG  = LOG Reflu2 ;

 K_gnie = A_gnie * (LOG (Reflu2 / 
               (( B_gnie * ReyLogG ) + C_gnie ))) ;
 F_frot = 0.5 / ( K_gnie ** 2. ) ;
 F_frot2 = F_frot ** 0.5 ;

 chpobor 'Reyn - 1000. ' ReyGnie ;
 chpobor 'coeff frott  ' ( 2. * F_frot );

 PrGnie = ((pranfil2 ** s_23 ) - 1. ) ;
 chpobor 'Pr**2/3 - 1  ' PrGnie ;
 nuss2 = F_frot * pranfil2 * ReyGnie / 
            ( 1. + (12.7 * F_frot2 * PrGnie )) ;
FINSI ;

SINON ;
 nuss1 = Nusslam ;
 nuss2 = Nusslam ;
FINSI ;

mess 'regime fluide 1 : ' RegEcoM.1 ;
mess 'regime fluide 2 : ' RegEcoM.2 ;

*** Coefficients d'echange
cech1  = nuss1 * lfil1 / ( TUBE1 . 'DH' )  ;
cech2  = nuss2 * lfil2 / ( TUBE2 . 'DH' )  ;


mess 'Diametres hydrauliques' ;
chpobor  'fluide 1     '  ( TUBE1 . 'DH' )   ;
chpobor  'fluide 2     '  ( TUBE2 . 'DH' )   ;

mess 'Conductivites ' ;
chpobor  'fluide 1  ' lfil1  ;
chpobor  'fluide 2  ' lfil2  ;

mess 'Nb de Prandtl ' ;
chpobor  'fluide 1  ' pranfil1;
chpobor  'fluide 2  ' pranfil2;

mess 'Nb de Nusselt ' ;
chpobor  'fluide 1  ' nuss1;
chpobor  'fluide 2  ' nuss2;

rvi . 'H1' =  'KCHT' $clim1 'SCAL' 'CENTRE' cech1 ;
rvi . 'H2' =  'KCHT' $clim2 'SCAL' 'CENTRE' cech2 ;
chpobor 'champ rvi echange 1  ' ( rvi . 'H1' ) ;
chpobor 'champ rvi echange 2  ' ( rvi . 'H2' ) ;

*** Proprietes physiques du solide
   tinter9 = 'NOEL' $VFW  (rvi . 'TSOL') ;
   lcb = ( pth_aci tinter9 ) . 'LAMS' ;
   rvi . 'Lacier' = 'KCHT' $VFW  'SCAL' 'CENTRE' lcb ;


MESS
'3bis) Activation des sources internes dans structures (en W/m3)' ;
***   3bis) Activation des sources internes dans les des structures (en W/m3)
*** Commentaires :
***  Comme lors de l'initialisation : on calcule la masse volumique de
***  chaque corps a partir de la temperature. Puis on en deduit les
***  puissances des sources internes.



MESS ' 4) Resolution de la thermique dans les structures ' ;
***   4) Resolution de la thermique dans les structures
   EXEC rvsol ;
opti echo 0 ;
 chpobor 'temperatures solides mini & maxi en Celsius (apres rvsol)'
        (( rvi.tsol ) - CK ) ;
*gounand opti echo 1 ;


MESS  '  5) Calcul du flux de chaleur aux interfaces Solides/Fluides' ;
***   5) Calcul du flux de chaleur aux interfaces Solides/Fluides

*** solide/Fluide 1
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $clim1 'OPER' 'MDIA' 'H1' 'INCO' 'TSOL' 'F1' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $clim1 'OPER' 'MDIA' 'H1' 'INCO' 'T1'  'F1' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   flx1 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                           ('NOMC' 'TSOL' (rvi . 'TSOL')))
                           ('KOPS' mmdia2 '*'
                           ('NOMC' 'T1'  (rvi . 'T1')))) ;
   'MESSAGE'
   'Flux integre sur interface solide/Fluide 1 (W) :'
            ('MAXIMUM' ('RESULT' flx1)) ;

*** solide/Fluide 2
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $clim2 'OPER' 'MDIA' 'H2' 'INCO' 'TSOL' 'F2' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $clim2 'OPER' 'MDIA' 'H2' 'INCO' 'T2'  'F2' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   flx2 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                           ('NOMC' 'TSOL' (rvi . 'TSOL')))
                           ('KOPS' mmdia2 '*'
                           ('NOMC' 'T2'  (rvi . 'T2')))) ;
   'MESSAGE'
   'Flux integre sur interface solide/Fluide 2 (W) :'
            ('MAXIMUM' ('RESULT' flx2)) ;


***   6) Transfert en une puissance volumique dans le fluide
   rvi . 'F1' = flx1 ; comm fluide un ;
   pflx1 = '*' mfc1 flx1 ;
   rvi . 'PFLU1' = 'KCHT' $mtflu1  'SCAL' 'CENTRE' pflx1 ;
  mess 'Puissance volumique fournie au fluide 1   : ' ;
  chpobor 'interface fluide 1 tube solide     ' pflx1 ;
chpobor ( rvi . 'PFLU1' ) ;
   rvi . 'F2' = flx2 ; comm fluide deux ;
   pflx2 = '*' mfc2 flx2 ;
   rvi . 'PFLU2' = 'KCHT' $mtflu2  'SCAL' 'CENTRE' pflx2 ;
  mess 'Puissance volumique fournie au fluide 2   : ' ;
  chpobor 'interface fluide 2 tube solide     ' pflx2 ;
chpobor ( rvi . 'PFLU2' ) ;

MESS '6bis) Activation des sources internes dans le fluide (en W/m3)' ;
***   6bis) Activation des sources internes dans le fluide (en W/m3)
*** Commentaires :
***  Comme lors de l'initialisation : on calcule la masse volumique de
***  chaque corps a partir de la temperature. Puis on en deduit les
***  puissances des sources internes.


MESS ' 7) Calcul des proprietes physiques fluide ' ;
***   7) Calcul des proprietes physiques fluide
***      a partir des variables a l'instant precedent.
***
*** Fluide 1 :
dmtflu1 = DOMA $mtflu1 'MAILLAGE' ;
Tflus = REDU dmtflu1 (rvi . 'TFLU') ;
Tfluc = 'NOEL' $mtflu1 (rvi . 'TFLU') ;

ppthe8s = PTH_HE8  Tflus ;
ppthe8c = PTH_HE8  Tfluc ;
rhes  = ppthe8s . 'ROG' ;
muhes = ppthe8s . 'MUG' ;
lhes  = ppthe8s . 'LAMG' ;
cphes = ppthe8s . 'CPG' ;
rhec  = ppthe8c . 'ROG' ;
muhec = ppthe8c . 'MUG' ;
lhec  = ppthe8c . 'LAMG' ;
cphec = ppthe8c . 'CPG' ;

uflu1s = ELNO $mtflu1 (KOPS ( TUBE1 . 'G' / Spass1) / rhec) ;
uflu1s = NOMC 'UZ' uflu1s ;
rvi . 'RCPFLU1' =  'KCHT' $mtflu1 'SCAL' 'CENTRE' ('*' rhec cphec) ;
rvi . 'LFLU1'   =  'KCHT' $mtflu1 'SCAL' 'CENTRE' lhec ;
rvi . 'UFLU1'   =  'KCHT' $mtflu1 'VECT' 'SOMMET' uflu1s ;

fb1 = (exco uz (rvi . 'UFLU1')) ;
chpobor 'rvi . uflu1   : ' fb1 ;

*** Fluide 2 :
dmtflu2 = DOMA $mtflu2 'MAILLAGE' ;
Tflus = REDU dmtflu2 (rvi . 'TFLU') ;
Tfluc = 'NOEL' $mtflu2 (rvi . 'TFLU') ;

ppthe8s = PTH_HE8  Tflus ;
ppthe8c = PTH_HE8  Tfluc ;
rhes  = ppthe8s . 'ROG' ;
muhes = ppthe8s . 'MUG' ;
lhes  = ppthe8s . 'LAMG' ;
cphes = ppthe8s . 'CPG' ;
rhec  = ppthe8c . 'ROG' ;
muhec = ppthe8c . 'MUG' ;
lhec  = ppthe8c . 'LAMG' ;
cphec = ppthe8c . 'CPG' ;

uflu2s = ELNO $mtflu2 (KOPS ( TUBE2 . 'G' / Spass2) / rhec) ;
uflu2s = NOMC 'UZ' uflu2s ;
rvi . 'RCPFLU2' =  'KCHT' $mtflu2 'SCAL' 'CENTRE' ('*' rhec cphec) ;
rvi . 'LFLU2'   =  'KCHT' $mtflu2 'SCAL' 'CENTRE' lhec ;
rvi . 'UFLU2'   =  'KCHT' $mtflu2 'VECT' 'SOMMET' uflu2s ;

fb2 = (exco uz (rvi . 'UFLU2')) ;
chpobor 'rvi . uflu2   : ' fb2 ;


MESS '  8) Calcul de la temperature dans l eau' ;
***   8) Calcul de la temperature dans l'eau
   EXEC rvflu ;
chpobor 'temperature fluide (Celsius) : '  ((rvi . 'TFLU') - CK );

*** Convergence :
   then = 'COPIER' (rvi . 'TFLU') ;
   tcbn = 'COPIER' (rvi . 'TSOL') ;
   errhe  = '/' ('MAXIMUM'  ('-' then thep) 'ABS')
                ('MAXIMUM'  thep 'ABS') ;
   errcb = '/' ('MAXIMUM'  ('-' tcbn tcbp) 'ABS')
               ('MAXIMUM'  tcbp 'ABS') ;
   'MESSAGE' ('CHAINE' 'Iteration ' &bouc) ;
   'MESSAGE' ('CHAINE' '   errhe  = ' errhe) ;
   'MESSAGE' ('CHAINE' '   errcb  = ' errcb) ;
   thep = then ;
   tcbp = tcbn ;
   testcvg = 'ET' (errhe '<' 1.D-4)  (errcb '<' 1.D-4) ;
   'SI' testcvg ;
      maxiter = faux ;
      'QUITTER' bouc ;
    'FINSI' ;

fin bouc ;

TCpu2 = tcpu9 0. ;


************************************************************************
*** **------------ FIN DU CALCUL  ------------------------------********
************************************************************************



***********************************************************************
***   POST TRAITEMENT SOMMAIRE
***********************************************************************


*** Calcul des Temperatures et les vitesses
*** ---------------------------------------
tk = rvi.tsol ; comm temperatures solides en K ;
tc = tk - ck ; comm temperatures solides en C ;
t_tn = rvsol . 'INCO' . 'TSOL' ;

TsoliC = (rvi . 'TSOL') - CK ;
TfluiC = (rvi . 'TFLU') - CK ;

*** Verification des bilans de puissance 
Wsolid = 'MAXIMUM' ('RESULT'(( rvi . 'Soursol' ) * 
        ( 'DOMA' $VFW 'VOLUME' ))) ; 

opti echo 0 ;

Ecform 'Diametre passage fluide 1           '
  '(F5.3)' 'm' ( TUBE1 . 'DH' ) ;
Ecform 'Diametre passage fluide 2           '
  '(F5.3)' 'm' ( TUBE2 . 'DH' ) ;

Ecform 'Surface de passage fluide 1         '
  '(E12.5)' 'm2' Spass1 ;
Ecform 'Surface de passage fluide 2         '
  '(E12.5)' 'm2' Spass2 ;

Ecform 'Debit massique fluide 1  '
  '(F5.2)' 'kg/s' ( TUBE1 . 'G' ) ;
Ecform 'Debit massique fluide 2  '
  '(F5.2)' 'kg/s' ( TUBE2 . 'G' ) ;


Ecform 'Volume solide          '
  '(E12.5)' 'm3' VoSol  ;

Ecform 'Densite de puissance interne          '
  '(E12.5)' 'W/m3' Sousol  ;
chpobor 'X_TOT (pour sources) : ' X_TOT ;

Ecform 'Puissance surfacique plasma          '
  '(E12.5)' 'W' W_SPLSM  ;
Ecform 'Puissance surfacique interne          '
  '(E12.5)' 'W' W_SBZ  ;
WsolTot = Wsolid + W_SPLSM + W_SBZ ;
WfluTot = ('MAXIMUM' ('RESULT' flx1)) + ('MAXIMUM' ('RESULT' flx2)) ;

Chpobor 'Coefficient echange fluide 1 (W/m2.K) ' Cech1 ;
Chpobor 'Coefficient echange fluide 2 (W/m2.K) ' Cech2 ;

ecform
'Flux interface solide/Fluide   1    '
 '(E12.5)'  'W'  ('MAXIMUM' ('RESULT' flx1)) ;
ecform
'Flux interface solide/Fluide   2    '
 '(E12.5)'  'W'  ('MAXIMUM' ('RESULT' flx2)) ;

ecform
'Puissance totale solide             '
 '(E12.5)'  'W'  WsolTot ;

ecform
'Puissance totale fluide             '
 '(E12.5)'  'W'  WfluTot ;

mess ;
chpobor 'temperatures solides min & max (Celsius)   '
  ((redu  ( rvi . 'TSOL' ) VFW) - CK ) ;
mess ;

chpobor 'temp. fluides min & max (Celsius) fluide 1     '
  ((redu  ( rvi . 'TFLU' ) tube1.mtflu ) - CK ) ;
chpobor 'temp. fluides min & max (Celsius) fluide 2     '
  ((redu  ( rvi . 'TFLU' ) tube2.mtflu ) - CK ) ;


si complet ;
 opti echo 1 ;
sinon ;
 opti echo 0 ;
finsi ;

si complet ;

*** Evolution fluide 1
*** ------------------
tfc9 = redu TfluiC tube1.mtflu ;
tfc9c = noel $mtflu1 tfc9 ;
mtfc9c = doma $mtflu1 centre ;
repe x (nbno mtfc9c )  ;
pt = mtfc9c point &x ;
val89 = coor 3 pt ;
wal89 = extr tfc9c scal pt ;
 si ( ega &x 1 ) ;
  dqqc = val89 ;
  drrc = wal89 ;
 sinon ;
  dqqc = dqqc et val89 ;
  drrc = drrc et wal89 ;
 finsi ;
fin x ;
mindqqc = mini dqqc ;
ndqqc = dime dqqc ;
dqqc = dqqc - ( prog ndqqc  *  mindqqc ) ;
evtfint1 =   evol manu  'Z' dqqc drrc  ;

*** Evolution fluide 2
*** ------------------
tfc9 = redu TfluiC tube2.mtflu ;
tfc9c = noel $mtflu2 tfc9 ;
mtfc9c = doma $mtflu2 centre ;
repe x (nbno mtfc9c )  ;
pt = mtfc9c point &x ;
val89 = coor 3 pt ;
wal89 = extr tfc9c scal pt ;
 si ( ega &x 1 ) ;
  dqqc = val89 ;
  drrc = wal89 ;
 sinon ;
  dqqc = dqqc et val89 ;
  drrc = drrc et wal89 ;
 finsi ;
fin x ;
mindqqc = mini dqqc ;
ndqqc = dime dqqc ;
dqqc = dqqc - ( prog ndqqc  *  mindqqc ) ;
evtfint2 =   evol manu  'Z' dqqc drrc  ;
*** Chercher les abscisses curvilignes

 trac cach TsoliC VFW titr 'Temperature solide' ;
 trac cach TsoliC sc1 titr 'Temperature interface solide-fluide 1' ;
 trac cach TsoliC sc2 titr 'Temperature interface solide-fluide 2' ;
 trac cach TfluiC tube1.mtflu titr 'Temperature fluide 1' ;
 trac cach TfluiC tube2.mtflu titr 'Temperature fluide 2' ;
 trac TfluiC ( tube1.mtflu et tube2.mtflu ) 
     titr 'Temperatures fluides 1 et 2 ' ;
 dess ( evtfint1 et evtfint2 ) mima
        titr 'Temperatures fluides 1 et 2 vs abscisse curviligne' ;
 dess EVQfw mima titr 'sources' ;

 factech = 2. ;

 loupe oeil TsoliC VFW 
     'Temperature solide' factech 0. ;
 loupe oeil TsoliC sc1 
     'Temperature interface solide-fluide 1' factech ;
 loupe oeil TsoliC sc2 
     'Temperature interface solide-fluide 2' factech ;
 loupe oeil TfluiC tube1.mtflu 
     'Temperature fluide 1' factech 5. ;
 loupe oeil TfluiC tube2.mtflu 
     'Temperature fluide 2' factech 5. ;
 loupe oeil TfluiC ( tube1.mtflu et tube2.mtflu ) 
     'Temperatures fluides 1 et 2' factech 0. ; 

finsi ;
************************************************************************
****************FIN DU POST-TRAITEMENT SOMMAIRE*************************
************************************************************************

'TEMPS' 'ZERO' ;

chpvide matvide = 'KOPS' 'MATRIK' ;
rvsol . METHINV.'MATASS'=matvide;
rvsol . METHINV.'MAPREC'=matvide;
rvflu . METHINV.'MATASS'=matvide;
rvflu . METHINV.'MAPREC'=matvide;
menage;
TCpu3 = tcpu9 0. ;

opti echo 0 ;
mess 'Duree calcul matrices des connectivites : ' TCpu1 ;
mess 'Duree calcul iteratif                   : ' TCpu2 ;
mess 'Duree menage                            : ' TCpu3 ;
opti echo 1 ;

FIN ;
 

© Cast3M 2003 - All rights reserved.
Disclaimer