Télécharger Th1D-T3D-Ebul.dgibi

Retour à la liste

* NOM DE FICHIER    : Th1D-T3D-Ebul.dgibi
 
*******************************************
*******************************************
 
************************************************************************
* DESCRIPTION       : Modelisation thermique 3D - thermohydraulique 1D avec possibilite
*                     d'ebullition nucleee sur la paroi chauffante en contact avec le
*                     fluide caloporteur (voir ci-dessous le domaine de calcul et les
*                     phenomenes physiques modelises).
*
**********************************************************************************************
*       A U T E U R S
**********************************************************************************************
*
* Julien BIGNALET   Stagiaire ENSEEIHT - Genie Atomique 2006 - mel : julien.bignalet@alumni.enseeiht.fr
* Jean-Louis LASZLO - CEA/DEN/DANS/DRSN/SIREN/LASPI - mel : jean-louis.laszlo@cea.fr
* Stephane GOUNAND  - CEA/DEN/DANS/DM2S/SFME/LTMF   - 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                            VRAI                  FAUX
*                  *******                            ****                  ****
*                                                                   Pas d'algorithme
*                                                                   "diphasique", donc
*                                                                   pas de seuillage
*                                                                   de la face externe
*                                                                   du crayon chauffant
*
* Pression                                            146.4 bar            146.4 bar
* Diametre externe crayon chauffant                     9.5 mm             9.5 mm
* Longueur chauffante                                 360.0 mm           360.0 mm
* Debit massique commun jambes chaude et froide       0.196 kg/s         0.196 kg/s
* Puissance lineique maximale crayon chauffant        480.0 W/cm         480.0 W/cm
* Flux critique (BOWRING) a mi-longueur chauffante    379.9 W/cm2        379.9 W/cm2
* Flux critique (BOWRING) a mi-longueur chauffante   1133.8 W/cm        1133.8 W/cm
* Puissance incidente crayon  chauffant             16094.0 W          16094.0 W
*
* Sources internes prises en compte
* Sources internes dans les structures
* Sources internes dans les fluides
* Sources internes structures                          15809.8 W          15791.8 W
* Sources internes fluides                               912.0 W            909.8 W
* Sources internes totales                             16721.8 W          16701.6 W
 
* Puissance interface
* jambe chaude & parois separant les 2 ecoulements       -370.4 W            581.1 W
* Puissance interface
* jambe froide & parois separant les 2 ecoulements       1589.8 W            638.4 W
* Puissance interface
* jambe froide & tubes de force                          2700.4 W           2083.3 W
* Pertes vers eau statique de la piscine                11600.8 W          13150.1 W
 
* Fluide caloporteur :
* ------------------
* Temperature entree fluide chaud (imposee)               252.0 degres C    252.0 degres C
* Temperature sortie fluide chaud                         269.5 degres C    270.6 degres C
* Gradient de temperature fluide chaud                     17.5 degres       18.6 degres
* Temperature entree fluide froid                         245.5 degres C    247.2 degres C
* Temperature sortie fluide froid (imposee)               251.0 degres C    251.0 degres C
* Gradient de temperature fluide froid                      5.5 degres        3.8 degres
 
* Temperature structures solides :
* ------------------------------
* Temperature de saturation                               340.2 degres C
* Temperature ebullition nucleee maximale (JENS & LOTTE)  342.9 degres C    PAS DE TRAITEMENT
* Temperature maximale paroi externe crayon chauffant     342.9 degres C    EBULLITION NUCLEEE
* Temperature maximale paroi externe en contact piscine    99.5 degres C     72.4 degres C
**********************************************************************************************
**********************************************************************************************
*
************************************************************************
* PARAMETRES PRINCIPAUX
************************************************************************
* discr    = |LINE : discretisation lineaire (suffisante)
*            |QUAF :  discretisation quadratique
* graph    =  VRAI : pour tracer des choses
* interact = |VRAI : pour une utilisation interactive avec traces a l'ecran
*            |FAUX : execution sans arret et sans traces
* rayo     = |VRAI : rayonnement (face a face) dans la cavite de gaz
*            |FAUX : non prise en compte du rayonnement
* verbose  =  FAUX : moins de messages dans le recapitulatif
 
* precis | = |VRAI : represente seulement les bonnes puissances
*        |   |FAUX : permet de comparer les differentes puissances aux
*        |           meme interfaces afin de choisir la plus precise
 
* egasol = |VRAI    : prise en compte sources internes reparties dans les structures
*          |FAUX    : pas de prise en compte sources internes reparties  dans les structures
* egaliq = |VRAI    : prise en compte sources internes reparties  dans les fluides
*          |FAUX    : pas de prise en compte sources internes reparties  dans les fluides
* varcna = |VRAI    : convection naturelle avec la piscine
*          |FAUX    : echange constant avec la piscine (valeur HCNCST)
 
* complet = |VRAI   : calcul avec toutes les iterations
*           |FAUX   : calcul rapide, avec nombre reduit d'iterations
************************************************************************
************************************************************************
*
*
***********************
* DOMAINE DE CALCUL : *
************************************************************************
*      PISCINE remplie d'eau statique               |--|----|------|---|
*                            --------              /                   |
*                                                 /                    |
* "mtef_t" ->|----------|----|--------|----|-----|                     |
*            |                                      |--|----|------|---|
*            |     tube force externe : "mt_tfe"   /                   |
*            |                                    /                    |
* "mten_t" ->|----------|----|--------|----|-----|                     |
*            |                                         |----|------|---|
*            |        lame de gaz (N2) : "mt_gaz"     /                |
*            |                                       /                 |
* "mtif_t" ->|----------|----|--------|----|-----|--|                  |
*            |                                                         |
*            |                                         |----|------|---|
*            |   tube de force interne : "mt_tfi"     /                |
*            |                                       /                 |
* "mtin_t" ->|----------|----|--------|----|-----|--|                  |
*            |                                        <==              |
*            |<-"bahe_p"                              <==    "hahe_p"->|
*            |            fluide froid : "mthe_p"     <==              |
*            |                                        <==              |
* "mtif_pe"->|----------|----|--------|----|-----|--|--|-----|-----|---|
*            |                  panier externe : "mt_pe"                       |
* "mtin_pe"->|----------|----|--------|----|-----|--|--|-----|-----|---|
*            |    lame d'eau stagnante : "mt_lam"                      |
* "mtif_p" ->|----------|----|--------|----|-----|--|--|-----|-----|---|
*            |                  panier interne : "mtcb_p"                      |
* "mtin_p" ->|----------|----|--------|----|-----|--|--|-----|-----|---|
*            |                                            ==>          |
*            |<-"bahe"    fluide chaud : "mthe"           ==>          |
*            |                                            ==>          |
* "mtif"   ->|----------|----|--------|----|-----|--|--|-----| "hahe"->|
*            |           gaine : "mtcb"                      |\        |
* "mtin"   ->|----------|----|--------|----|-----|--|--|-----| \       |
*          ->| rallonge |bou | zone   |bou | rallonge sup :  |  \      |
*            |   inf :  |chon| chauf. |chon|                 |   \     |
*            | "mtcbAr" |    |        |    |    "mtcbEr"     |    \    |
*   Axe    --------------------------------------------------------|---|->
* symetrie -------------------------------------------------------------->
*            |          |    |        |    |     |  |  |     |     |   |
*
* Parties :  |   A      | B  |   C    |  D |  E  |F |G | H,I |  J  | K |
*            |          |    |        |    |     |  |  |     |     |   |
*            z1        z2   z3       z4   z5    z6 z7 z8  z9 z10  z11  zM
*
*
*--------------------------------|
* PHENOMENES PHYSIQUES MODELISES |
*--------------------------------|
*
*
*                  PISCINE : Temperature fixe       |------------------|
*                                                  /                   |
*                                                 /                    |
*            |--------Convection naturelle-------|                     |
*            |                                      |------------------|
*            | tube force externe : Conduction     /                   |
*            |                                    /                    |
*            |-----------------------------------|                     |
*            |                                         |---------------|
*            | lame de gaz (N2) :                     /                |
*            |          Rayonnement + Conduction     /                 |
*            |--------------------------------------|                  |
*            |                                                         |
*            |                                         |----|------|---|
*            | tube de force interne : Conduction     /                |
*            |                                       /                 |
*            |---- Convection forcee ---------------|                  |
*            |  <==                                                    |
*   Temp     |  <==                                                    |
*   imposee  |  <==          fluide froid                              |
*            |  <==                                                    |
*            |---- Convection forcee ----------------------------------|
*            |                 panier externe : Conduction             |
*            |---------------------------------------------------------|
*            |                     lame d'eau : Conduction             |
*            |---------------------------------------------------------|
*            |                 panier interne : Conduction             |
*            |---- Convection forcee ----------------------------------|
*   Temp     |                                            ==>          |
*   imposee  |              fluide chaud                  ==>          |
*            |                                            ==>          |
*            |---- Convection forcee ------------------------|         |
*            |                          gaine : Conduction   |\        |
*            |-----------------------------------------------| \       |
*            | rallonge |bou | crayon |bou | rallonge sup :  |  \      |
*            |   inf :  |chon|        |chon|                 |   \     |
*            |Conduction|vide| chauff.|vide|   Conduction    |    \    |
*   Axe    --------------------------------------------------------|---|->
* symetrie -------------------------------------------------------------->
 
************************************************************************
*****************DECLARATIONS DES PROCEDURES****************************
************************************************************************
 
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 ;
 
debproc recal coul0*mot dhp0*chpoint mod0/MMODEL m0*maillage
            pt1*point pt2*point npt1*entier
           pt3/point npt2/entier ;
si ( non ( existe mod0)) ;
 chp0 =  dhp0 ;
sinon ;
 chp0 = ELNO mod0 dhp0 ;
 sele0 = DOMA mod0 'XXVOLUM' ; comm surface elements ;
 densf0 = KOPS '/' dhp0 sele0 ; comm chpo scal centre ;
 densf1 = ELNO mod0 densf0 ; comm chpo scal sommet ;
finsi ;
qt1 = m0 poin proc pt1 ;
qt2 = m0 poin proc pt2 ;
opti elem seg2 ;
si (existe pt3 );
 qt3 = m0 poin proc pt3 ;
 drt0 = droi  npt1 qt1 ( droi npt2 qt2 qt3 ) ;
sinon ;
 drt0 = droi npt1 qt1 qt2 ;
finsi ;
* * castem13 et avant
* si ( non ( existe mod0)) ;
*  tdrt0 = INT_COMP m0 chp0 drt0 ;
* sinon ;
*  tdrt0 = INT_COMP m0 densf1 drt0 ;
* finsi ;
***** modif bp du 03/09/2013 suite a l'evolution #7817 *****
si ( non ( existe mod0)) ;
 chp0ok = chp0 ;
sinon ;
 chp0ok = densf1;
finsi ;
* normale aux elements (d'amplitude=0.1*taille des elements)
modcoq4 = mode m0 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'COQ4';
vcoq4 = VSUR modcoq4;
vcoq4x = exco vcoq4 'VX' 'SCAL';
vcoq4y = exco vcoq4 'VY' 'SCAL';
vcoq4z = exco vcoq4 'VZ' 'SCAL';
nvcoq4 = ((vcoq4x**2) + (vcoq4y**2) + (vcoq4z**2))**0.5;
vcoq4 = 0.1 * (nvcoq4**-0.5) * 
       (vcoq4 exco (mots VX VY VZ) (mots UX UY UZ));
*trac (vect 1. modcoq4 vcoq4 (mots UX UY UZ)) m0;
vcoq4 = chan vcoq4 'CHPO' modcoq4 'MOYENNE';
* extrusion pour creation elements massifs 
* (seuls acceptes par PROI qui est appele par INT_COMP)
m0pp chp0pp = m0 chp0ok plus vcoq4;
m0mm chp0mm = m0 chp0ok moin vcoq4;
opti elem CUB8;
m0ok = (volu m0mm m0) et (volu m0 m0pp);
chp0ok = chp0mm et chp0ok et chp0pp;
*trac chp0ok m0ok;
tdrt0 = INT_COMP m0ok chp0ok drt0 ;
***** fin modif bp du 03/09/2013 suite a l'evolution #7817 *****
ev0 = evol coul0 'CHPO' tdrt0 drt0 ;
si ( non ( existe mod0)) ;
 respro ev0 ;
sinon ;
 respro ev0 densf1 ;
finsi ;
finproc ;
 
DEBPROC BOWRING TAB_1*TABLE chlat/flottant ;
*
*
* CALCUL DU FLUX CRITIQUE SUIVANT LA CORRELATION DE BOWRING
* --- entrees
*
INIVEAU1 = TAB_1.'NIVEAU' ;
D_DIAM1  = TAB_1.'D_DIAM' ;
L_LONG1  = TAB_1.'L_HEATED' ;
P_PRES1  = TAB_1.'P_IN' ;
G_VITE1  = TAB_1.'G_IN' ;
T_TEMP1  = TAB_1.'T_IN'  ;
 
 
* - test sur la vitesse de l'eau
SI ((G_VITE1  < 136.) OU ( G_VITE1 > 18600.))                    ;
 MESS 'Vitesse massique : ' G_VITE1;
 ERRE '@BOWRING -> Vitesse massique hors [136. , 18600.] (Kg/M2/S)';
FINSI                                                                 ;
 
* - test sur le diametre
SI ((D_DIAM1 < 2.E-3) OU (D_DIAM1 > 450.E-3))                       ;
 MESS 'Diametre : '  D_DIAM1;
 ERRE '@BOWRING -> Diametre hors [0.002 0.45] (M)'             ;
FINSI                                                                 ;
 
* - test sur la Pression
SI ((P_PRES1 < 1.E5) OU (P_PRES1 > 200.E5))                      ;
 MESS 'Pression : ' P_PRES1;
 ERRE '@BOWRING -> Pression hors de [1.E5, 200.E5] (Pa) '       ;
FINSI                                                                 ;
 
* - test sur la longueur chauffee
SI ((L_LONG1 < 0.15) OU (L_LONG1 > 3.7))                    ;
 MESS 'Longueur : ' L_LONG1;
 ERRE '@BOWRING --> Longueur hors de [0.15,3.7](M) '  ;
FINSI                                                                 ;
 
P1 = P_PRES1 / 6900000. ;
 
SI (P1 &lt;EG 1.) ;
 F1 = (((P1 ** 18.942) * (EXP (20.8 * (1. - P1)))) + 0.917) / 1.917;
 F2 = (F1 * 1.309)/(((P1 ** 1.316)*(EXP(2.444*(1. - P1)))) + 0.309);
 F3 = (((P1 ** 17.023)*(EXP(16.658*(1. - P1)))) + 0.667)/1.667;
 F4 = F3 * (P1 ** 1.649) ;
SINON ;
 F1 = (P1 ** (-0.368))*(EXP(0.648*(1. - P1)));
 F2 = (P1 ** (-0.448))*(EXP(0.245*(1. - P1)));
 F3 = P1 ** 0.219;
 F4 = F3 * (P1 ** 1.649) ;
FINSI ;
 
 
si (existe chlat) ;
 L_VAP = chlat ;
sinon ;
 L_VAP = 'VARI' 'LATENT' P_PRES1 ;
finsi ;
 
S_SAT = hlsat - ( 'VARI' 'HLS' P_PRES1 ( T_TEMP1 + ck ) ) ;
 
A__1 = 0.5793 * L_VAP * D_DIAM1 * G_VITE1 * F1 /
     (1. + (0.0143 * F2 * (D_DIAM1 ** .5) * G_VITE1 )) ;
 
B__1 = .25 * D_DIAM1 * G_VITE1 ;
 
C__1 = 0.077 * D_DIAM1 * G_VITE1 * F3 /
     (1. + (0.347 * F4 * ((G_VITE1/1356.) ** (2. - (.5 * P1))))) ;
 
QCHFW = (A__1 + (B__1 * S_SAT)) / (C__1 + L_LONG1) ;
 
G1 = G_VITE1 * PI * D_DIAM1 * D_DIAM1 / 4. ;
*
FINPROC   QCHFW   ;
 
 
 
debproc loupe m0*maillage c0*chpoint gross0*flottant
    tjt0*mot tipiso*mot cnt0/entier ;
opti isov tipiso ;
eye = 1.e30 1.e30 1.e30 ;
* 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 (existe cnt0) ;
 si ( cnt0 < 0 );
  trac cach eye c0 m0 (cont m0 )
  coupe (0. 0. 0.) (0. 1. 0.) (0. 0. 1.) titr tjt0 ;
 sinon ;
*  trac cach eye c0 m0 (cont m0 ) titr tjt0 ;
  trac cach eye c0 m0 titr tjt0 ;
 finsi ;
sinon ;
trac cach eye c0 m0 titr tjt0 ;
finsi ;
'FORME' orig0 ;
finproc ;
 
 
DEBPROC h2o_psat t*'FLOTTANT' ;
* H2O_PSAT  Saturation pressure at a given temperature in kPa
* H2O_PSAT(T) Returns the saturation pressure at a given
* temperature.  Based on the correlation (K function) given in the
* ASME STEAM TABLES - SIXTH EDITION
* Temperature range from 0.0�C to 374.15�C
*  Called function: NONE
*  Required Inputs are: T - temperature in �C
* ---------------------------------------------------------------
* The MATLAB function was created by Tibor Balint, December 1998
* TBoreal Research Corporation, Toronto, Ont. Canada
* (tibor@netcom.ca) and also, University of Warwick, UK
* ---------------------------------------------------------------
 
 
SI  ( T > 374.15 ) ;
 MESS 'Temperature out of range (above 374.15�C)' ;
FINSI ;
SI ( T < 0.0 ) ;
 MESS 'Temperature out of range (below 0.0�C)' ;
FINSI ;
 
 
T1=T+273.15;
* set up the constants
* the reduced temperature is calculated from theta=T/Tc1
* temperature constants used in the reduced temperature Tc1=647.3K
TC1=647.3;
*the reduced saturation pressure is calculated from
* betaK(theta)=ps/pc1, where ps=ps(T)
*pc1=22120000 N/m^2
pc1=22120000;
 
*Saturation line constants from k1 to k9
k = TABLE ;
s = TABLE ;
k . 1 = -7.691234564e0 ;
k . 2 =    -2.608023696e1 ;
k . 3 =    -1.681706546e2 ;
k . 4 =    6.423285504e1 ;
k . 5 =    -1.189646225e2 ;
k . 6 =    4.167117320e0 ;
k . 7 =    2.097506760e1 ;
k . 8 =    ( 10 ** 9. ) ;
k . 9 =    6. ;
 
* build up the correlation
a1 = 0. ;
 
REPE v 5 ;
   s . &v = ( k . &v ) * ((1.-(T1/TC1)) ** &v ) ;
   a1 = a1 + (s . &v) ;
FIN v  ;
 
a2 = 1. / (T1/TC1) ;
a3 = (1. - (T1/TC1));
a4 = 1. + ((k . 6) * a3) + ((k . 7) * (a3 ** 2.)) ;
a5 = (k . 8) * (a3 ** 2.) + (k . 9) ;
 
*the reduced saturation pressure is
betaK = exp(a2*(a1/a4)-(a3/a5));
 
PSAT = betaK * pc1 ; comm saturation pressure in Pa ;
 
FINPROC psat ;
 
DEBPROC h2o_tsat p*'FLOTTANT' ;
* H2O_TSAT  Saturation temperature at a given pressure in �C
* H2O_TSAT(P) Returns the saturation temperature at a given
* pressure.  Based on the correlation (K function) given in the
* ASME STEAM TABLES - SIXTH EDITION
* Pressure range from 0.6108kPa to 22120kPa
*  Called function: h2o_psat
*  Required Inputs are: P - pressure in Pa
* ---------------------------------------------------------------
*  The MATLAB function was created by Tibor Balint, December 1998
* TBoreal Research Corporation, Toronto, Ont. Canada
* (tibor@netcom.ca) and also, University of Warwick, UK
* ---------------------------------------------------------------
 
 
SI ( P > 22120.e3 ) ;
 MESS 'Pressure out of range (above 22120 kPa)' ;
FINSI ;
SI ( P < 0.61080062637844 )  ;
 MESS 'Pressure out of range (below 0.6108 kPa)' ;
FINSI;
 
PIN=P;
*set the minimum temperature to the lower range of the h2o_psat function
TMIN=0;
*set the maximum temperature to the lower range of the h2o_psat function
TMAX=374.15;
tm = 0.5 * ( tmin + tmax ) ;
pm =( h2o_psat tm);
 
* Iterate for the pressure by guessing the temperature and
* compare the result to the requested pressure input
REPE fb ;
 
SI (( abs (PIN - Pm)) > 0.1 ) ;
*SI ( &fb < 10 ) ;
 
   SI ( PIN &lt;EG pm ) ;
      TMAX = tm ;
 
   SINON ;
      TMIN = tm ;
   FINSI  ;
   tm = 0.5 * ( tmin + tmax ) ;
 
   Pm = ( h2o_psat tm)  ;
 
SINON ;
 QUITTER fb ;
FINSI ;
 
FIN fb ;
 
* Assign the last iteration value to the saturation temperature
TSAT=Tm ;
 
FINPROC tsat ;
 
DEBPROC h2o_cplq t*'FLOTTANT' p*'FLOTTANT' ;
 
TS = ( h2o_tsat p ) ; comm  calculate the saturation temperature ;
PS = ( h2o_psat T ) ; comm calculate the saturation pressure ;
 
SI ( T > TS ) ;
    MESS 'Cp : Temperature is above saturation for the given pressure' ;
FINSI ;
SI ( PS > P ) ;
   MESS 'Cp : Pressure is above saturation for the given temperature' ;
FINSI ;
 
 hl2 = ((((0.3654e-12 * p ) - 0.30413e-5 ) * p ) + 0.40047e4 ) ;
 hl3 = ((((-0.95261e-8 * p ) - 0.25785 ) * p ) + 0.20641e8 ) ;
 hl4 = 399.98;
 CPLIQ = hl2 + ( hl3 / ((hl4 - t ) ** 2.0)) ;
FINPROC cpliq ;
 
DEBPROC h2o_mu9 t*'FLOTTANT' p*'FLOTTANT' ;
 
* H2O_MU  Dynamic viscosity at a given temperature and density
* H2O_MU(T,RHO) Returns the dynamic viscosity at a given temperature and density.
* Dynamic viscosity in (Pa s) i.e., (kg/ms)
* Based on the correlation given in Appendix 6 of the
* ASME STEAM TABLES - SIXTH EDITION
* Range of validity of equation:
*    P<=500MPa for   0�C <= T <= 150�C
*    P<=350MPa for 150�C <= T <= 600�C
*    P<=300MPa for 600�C <= T <= 900�C
* Called function: none
*  Required Inputs are: T   - temperature in �C
*
*
* ---------------------------------------------------------------
* The MATLAB function was created by Tibor Balint, December 1998
* TBoreal Research Corporation, Toronto, Ont. Canada
* (tibor@netcom.ca) and also, University of Warwick, UK
* ---------------------------------------------------------------
 
*check the pressure and temperature ranges
TS = ( h2o_tsat p ) ; comm  calculate the saturation temperature ;
PS = ( h2o_psat T ) ; comm calculate the saturation pressure ;
 
 
SI ( T > TS ) ;
    MESS 'MU : Temperature is above saturation for the given pressure' ;
FINSI ;
SI ( PS > P ) ;
   MESS 'MU : Pressure is above saturation for the given temperature' ;
FINSI ;
 
* Get the saturated fluid density
RHO = ( h2o_rhol T P ) ;
 
 
*Reference temperature 647.27K
Tstar=647.27;
*Reference density 317.763 kg/m^3
RHOstar=317.763;
*Reference pressure 22.115x10^6 Pa
Pstar=22.115e6;
*Reference viscosity 55.071x10^-6 Pa s (or kg/ms)
mustar=55.071e-6;
 
tab = TABLE ;
tab . 1  = TABLE ;
tab . 2  = TABLE ;
tab . 3  = TABLE ;
tab . 4  = TABLE ;
tab . 5  = TABLE ;
tab . 6  = TABLE ;
tab . 7  = TABLE ;
 
 
*Coefficients H_i for  mu_0
HI = TABLE ;
HI . 1 = 1. ;
HI . 2 = 0.978197 ;
HI . 3 = 0.579829 ;
HI . 4 = -0.202354 ;
*list HI ;
*HI=[1 0.978197 0.579829 -0.202354];
 
*Coefficients H_ij for  mu_1
*The coefficients not listed are all 0, so first set the matrix elements to 0
REPE i 6 ;
  REPE j 7 ;
  tab . &i . &j = 0. ;
  FIN j ;
FIN i ;
*list tab ;
*list ( tab . 1 ) ;
*HIJ(1:6,1:7)=0;
 
H00=0.5132047; comm    ;
tab . 1 . 1 = H00;
H10=0.3205656; comm    ;
tab . 2 . 1 = H10;
 
 
 
H40=-0.7782567; comm    ;
tab . 5 . 1 =H40;
H50=0.1885447; comm    ;
tab . 6 . 1 =H50;
H01=0.2151778; comm    ;
tab . 1 . 2 =H01;
H11=0.7317883; comm    ;
tab . 2 . 2 =H11;
H21=1.241044; comm    ;
tab . 3 . 2 =H21;
H31=1.476783; comm    ;
tab . 4 . 2 =H31;
H02=-0.2818107; comm    ;
tab . 1 . 3 =H02;
H12=-1.070786; comm    ;
tab . 2 . 3 =H12;
H22=-1.263184; comm    ;
tab . 3 . 3 =H22;
H03=0.1778064; comm    ;
tab . 1 . 4 =H03;
H13=0.4605040; comm    ;
tab . 2 . 4 =H13;
H23=0.2340379; comm    ;
tab . 3 . 4 =H23;
H33=-0.4924179; comm    ;
tab . 4 . 4 =H33;
H04=-0.04176610; comm    ;
tab . 1 . 5 =H04;
H34=0.1600435; comm    ;
tab . 4 . 5 =H34;
H15=-0.01578386; comm    ;
tab . 2 . 6 =H15;
H36=-0.003629481; comm    ;
tab . 4 . 7 =H36;
*list tab ;
 
T1 = T + 273.15 ; comm    Celsius to K ;
Tbar = T1 / Tstar ;
* In the viscosity equation the firs (mu0) term of the product gives the viscosity
* of the steam in the ideal-gas limit and calculated from:
comp1=0;
REPE i 4 ;
   comp1 = comp1 + ( ( hi . &i  ) / (Tbar ** (&i-1)) );
FIN i ;
mu0 = (((Tbar) ** 0.5 ) / comp1 ) ;
* The second multiplicative factor is
RHObar = RHO / RHOstar ;
comp2=0;
REPE i 6 ;
REPE j 7 ;
 
      comp2 = comp2 + ( ( tab . &i . &j )
      * ( ((1. / Tbar ) - 1. ) ** ( &i - 1. ) )
      * ((RHObar - 1.) ** ( &j - 1. )) ) ;
FIN  j ;
FIN i ;
 
mu1 = exp ( RHObar * comp2 );
 
* The reduced viscosity is
mubar = mu0 * mu1 ;
 
* The dynamic viscosity is calculated from the reduced and reference viscosities
MU = mubar * mustar ;
 
 
 
FINPRO mu ;
 
DEBPROC h2o_lmda t*'FLOTTANT' p*'FLOTTANT' ;
hl = VARI 'HLS' p ( t + 273.15 ) ;
ak0=5.73738622e-1;
ak1=2.536103551e-1;
ak2=-1.45468269e-1;
ak3=1.387472485e-2;
xk = ( hl*1.7196904e-6 ) ;
lmda  = (((((ak3*xk)+ak2)*xk)+ak1)*xk)+ak0;
 
 
FINPROC lmda ;
 
 
DEBPROC h2o_rhol ;
ARGU t*'FLOTTANT' ; comm temperature in �C ;
ARGU p*'FLOTTANT' ; comm pressure in Pa ;
 
* ---------------------------------------------------------------
* H2O_RHOL  Subcooled density at a given temperature and pressure
* H2O_RHOL(T,P) Returns the subcooled density a given temperature
* and pressure.  Based on the correlation (Subsection 1) given in the
* ASME STEAM TABLES - SIXTH EDITION
* Temperature range from 0.0�C to 374.15�C
* Pressure range from 0.6108kPa to 22120kPa
*  Called function: h2o_psat(T), h2o_tsat(P)
*  The MATLAB function was created by Tibor Balint, December 1998
* TBoreal Research Corporation, Toronto, Ont. Canada
* (tibor@netcom.ca) and also, University of Warwick, UK
* ---------------------------------------------------------------
 
 
*check the pressure and temperature ranges
TS=h2o_tsat P ; comm calculate the saturation temperature ;
PS=h2o_psat T ; comm calculate the saturation pressure ;
 
si ( T > TS ) ;
   mess 'RHOL : Temperature is above saturation for the given pressure' ;
SINON;
 mess'RHOL : Pressure is above saturation for the given temperature' ;
finsi;
 
T1=T+273.15; comm convert the temperature from �C to K ;
 
* set up the constants
* the reduced temperature is calculated from theta=T/Tc1
* temperature constants used in the reduced temperature Tc1=647.3K
TC1=647.3; comm  the reduced temperature ;
theta=T1/TC1;
 
* The reduced pressure is: beta=ps/pc1
*where:  pc1=22120000 N/m^2
pc1=22120000.D0; comm reduced pressure i;
beta=P/pc1;
 
* Primary constants for Subregion 1
A11=7.982692717e0;
A12=-2.616571843e-2;
A13=1.522411790e-3;
A14=2.284279054e-2;
A15=2.421647003e2;
A16=1.269716088e-10;
A17=2.074838328e-7;
A18=2.174020350e-8;
A19=1.105710498e-9;
A20=1.293441934e1;
A21=1.308119072e-5;
A22=6.047626338e-14;
a1=8.438375405e-1;
a2=5.362162162e-4;
a3=1.72e0;
a4=7.342278489e-2;
a5=4.975858870e-2;
a6=6.5371543e-1;
a7=1.15e-6;
a8=1.1508e-5;
a9=1.4188e-1;
ab10=7.002753165e0;
ab11=2.995284926e-4;
ab12=2.04e-1;
 
Y=1 - (a1*(theta**2.)) - (a2*(theta**(-6.)));
Z=Y + (((a3*(Y**2.)) - (2.*a4*theta) + (2.*a5*beta))**0.5) ;
 
f1 = (A11*a5*(Z**(-5./17.))) ;
f2 = ( A12 + (A13*theta) + (A14*(theta**2.)) + (A15*((a6-theta)**10.))
+ ( A16*((a7 + (theta**19.))**(-1.)) ) );
f3 = (-1.) * ((a8 + (theta**11.))**(-1.))*(A17 + (2.*A18*beta) +
    (3.*A19*(beta**2.))) ;
f4 = (-1. ) *  (A20*(theta**18.)*(a9 + (theta**2.))*
        ((-3.*((ab10+beta)**(-4.))) + ab11)) ;
f5 =  (3.*A21*(ab12-theta)*(beta**2.));
f6 =  (4.*A22*(theta**(-20.))*(beta**3.)) ;
* The reduced volume is
x1 = f1 + f2 + f3 + f4 + f5 + f6 ;
 
* The reduced volume for subcooled fluid in Subregin 1 is
* calculated from x1=v/vc1 where vc1=0.00317 m^3/kg
vc1=0.00317;
vl=x1*vc1;
RHOL=1/vl;
FINPROC rhol ;
 
 
 
DEBPROC PRPN2 TS/'FLOTTANT' T/'CHPOINT' ;
 
***  Description
***  ***********
 
* Procedure donnant les proprietes du diazote dans la cavite de gaz :
*    - la masse volumique en approximant N2 a un gaz parfait
*    - la conductivite a partir des tables des Techniques de l'Ingenieur
* Loi des GP : PV = nRT donne P/rho = rT <=> rho = P/rT
* On est a une presssion de 5 bars, soit 0.5 MPa.
* P_N2 : pression de N2 en bars
* P   : pression en Pa
* R1  : constante des GP
* M   : masse molaire de N2 en kg/mol
* rho : masse volumique en kg/m3
* pres : pression de la table en MPa
* temp : temperature de la table en K
* ld_i : conductivite thermique de 200 a 1000 K a la pression i,
*        en mW/m/K
*
* T : temperature, parametre d'entree en Kelvin
*
 
 
P_N2 = 5. ;
P_MPa = P_N2 / 10. ;
 
**Masse volumique AZOTE
P = P_N2 * 1.e5 ;
R1 = 8.314 ;
M = 28.0134 * 1.e-3 ;
r = R1/M ;
 
**Conductivite AZOTE
pres  = prog 0.1 1. 2. ;
temp  = prog 200. 250. 300. 350. 400. 450. 500. 600. 700. 800. 1000. ;
ld_01 = prog 18.6 22.5 26.1 29.4 32.7 35.8 38.9 44.8 50.7 56.4 67.6 ;
ld_1  = prog 19.2 23.4 26.9 30.1 33.3 36.4 39.4 45.3 51.1 56.8 67.9 ;
ld_2  = prog 19.9 23.4 26.9 30.1 33.3 36.4 39.4 45.3 51.1 56.8 67.9 ;
 
nb1 = 11 ;
ld_05 = prog ;
 
REPETER bouc1 nb1 ;
  ldb_01 = EXTR ld_01 &bouc1 ;
  ldb_1  = EXTR ld_1 &bouc1 ;
  ldb_2  = EXTR ld_2 &bouc1 ;
  ld_T   = prog ldb_01 ldb_1 ldb_2 ;
  ldb_05 = IPOL P_MPa pres ld_T ;
  ld_05  = ld_05 ET (prog ldb_05) ;
FIN bouc1 ;
 
ld_05 = ld_05 * 1e-3 ;
 
**Resultats :
SI (EXISTE TS) ;
  rho = P / (r*TS) ;
  ld_CH = IPOL TS temp ld_05 ;
FINSI ;
 
SI (EXISTE T) ;
  rho = KOPS P / (KOPS r * T) ;
  ld_CH = IPOL T temp ld_05 ;
FINSI ;
 
RESPROC rho ld_CH ;
 
FINPROC ;
 
 
 
************************************************************************
* NOM         : PRPEAU
*-----------------------------------------------------------------------
*                Proprietes thermiques de l'eau a 150 bars
*-----------------------------------------------------------------------
* Input : PTOT : pression en Pa ou en bars
*         T ou TS : temperature sous forme de CHPOINT ou LISTREEL ou
*         FLOTTANT (en K)
*
* Output : RHOEAU : densite de l helium (USI) de meme type que T ou TS
*          MUEAU  : viscosite
*          LEAU   : conductivite thermique
*          CPEAU  : capacite calorifique a pression constante
*-----------------------------------------------------------------------
'DEBPROC' PRPEAU prt7*'LISTREEL' roo*'LISTREEL' muo*'LISTREEL'
        cpo*'LISTREEL' lambo*'LISTREEL' T/'CHPOINT'
        TL/'LISTREEL'  TS/'FLOTTANT' ;
 
   'SI' ('EXISTE' T) ;
      rhoeau = 'IPOL' T prt7 roo ;
      mueau  = 'IPOL' T prt7 muo ;
      leau   = 'IPOL' T prt7 lambo ;
      cpeau  = 'IPOL' T prt7 cpo ;
      'RESPRO' rhoeau mueau leau cpeau ;
   'FINSI' ;
   'SI' ('EXISTE' TL) ;
      rhoeau = 'IPOL' TL prt7 roo ;
      mueau  = 'IPOL' TL prt7 muo ;
      leau   = 'IPOL' TL prt7 lambo ;
      cpeau  = 'IPOL' TL prt7 cpo ;
      'RESPRO' rhoeau mueau leau cpeau ;
   'FINSI' ;
   'SI' ('EXISTE' TS) ;
      rhoeau = 'IPOL' TS prt7 roo ;
      mueau  = 'IPOL' TS prt7 muo ;
      leau   = 'IPOL' TS prt7 lambo ;
      cpeau  = 'IPOL' TS prt7 cpo ;
      'RESPRO' rhoeau mueau leau cpeau ;
   'FINSI' ;
*tc mise en commentaire du finsi ci dessous
*'FINSI';
*
* End of procedure file PRPEAU
*
'FINPROC' ;
 
 
 
'DEBPROC' PRPZIRC T/'CHPOINT' TL/'LISTREEL' TS/'FLOTTANT' ;
 
* Permet de donner la conductivite thermique du Zircaloy 4
* et la masse volumique.
* rho20 : masse volumique du zirconium a 20 degres C
* dil20 : coefficient de dilatation lineique
* Domaine de validite :| 0 a 1200 degres C pour la conductivite
*                      | 20 a 700 degres C pour la masse volumique
*
*      X.2. Arguments d'entree et de sortie
* Input = T : temperature sous forme de CHPOINT ou FLOTTANT (en K)
* Output = LZIRC : conductivite thermique du Zircalloy (W.m-1.K-1)
*          RHOZIRC : masse volumique (kg/m3)
 
*      X.3. Ecriture ==> voir isa4reb.procedur
 
***Pour la conductivite
prt = 'PROG' 273.15 'PAS' 100. 'NPAS' 12 ;
prk_zy4 = PROG 12.650 13.478 14.443 15.545 16.784 18.160
               19.673 21.324 23.111 25.036 27.098 29.297 31.633 ;
***Pour la masse volumique
rho20 = 6550. ;
dil20 = 6.e-6 ;
T20 = 293. ;
 
***
'SI' ('EXISTE' T) ;
   lzirc = 'IPOL' T prt prk_zy4 ;
   DT = T - T20 ;
   rhozirc = KOPS rho20 / ((1. + (DT * dil20)) ** 3.) ;
   'RESPRO' rhozirc lzirc ;
'FINSI' ;
'SI' ('EXISTE' TL) ;
   lzirc = 'IPOL' TL prt prk_zy4 ;
   DT = TL - T20 ;
   rhozirc = rho20 / ((1. + (DT * dil20)) ** 3.) ;
   'RESPRO' rhozirc lzirc ;
'FINSI' ;
'SI' ('EXISTE' TS) ;
   lzirc = 'IPOL' TS prt prk_zy4 ;
   DT = TS - T20 ;
   rhozirc = rho20 / ((1. + (DT * dil20)) ** 3.) ;
   'RESPRO' rhozirc lzirc ;
'FINSI' ;
*'FINSI' ;
'FINPROC' ;
 
 
 
************************************************************************
* NOM         : TRACCHPO
* DESCRIPTION : Trace d'un chpoint avec titre optionnel.
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stephane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mel : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 14/10/2002, version initiale
* HISTORIQUE : v1, 14/10/2002, creation
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Priere de PRENDRE LE TEMPS de completer les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'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' ('NON' ('EXISTE' mt2)) ;
   mt2 = mt ;
'FINSI' ;
'SI' ('EXISTE' oeil) ;
      'TRACER' 'CACH' oeil rescal mt mt2 'TITR' titpn nbisov ;
'SINON' ;
      'TRACER' 'CACH' rescal mt mt2 'TITR' titpn nbisov ;
'FINSI' ;
*
* End of procedure file TRACCHPO
*
'FINPROC' ;
 
 
debproc rayon0 m*maillage ;
x = coor 1 m ;
y = coor 2 m ;
r = ((x * x ) + ( y * y )) ** 0.5 ;
finproc r ;
 
 
DEBPROC VOLUJU sfluc*'MAILLAGE' sfluf*'MAILLAGE'
               licb*'MAILLAGE' lecb*'MAILLAGE' scb*'MAILLAGE'
               lipa*'MAILLAGE' lepa*'MAILLAGE' spa*'MAILLAGE'
               lipe*'MAILLAGE' lepe*'MAILLAGE' spe*'MAILLAGE'
               litf*'MAILLAGE' letf*'MAILLAGE' stf*'MAILLAGE'
               zn*'FLOTTANT' zn1*'FLOTTANT' nmail*'ENTIER'
               tabcb*'TABLE' tabpa*'TABLE' tabpe*'TABLE' tabtf*'TABLE'
               tabfluc*'TABLE' tabfluf*'TABLE' itab*'ENTIER' ;
 
** Description :
** ***********
* Procedure permettant de construire les maillages (volume) de chaque
* partie axiale du dispositif ainsi que les tables de connectivite
* entre fluide et solides.
* Les volumes 3D sont extrudees a partir des surfaces et des lignes 2D
* perpendiculaires a l'axe z.
*
* sfluc : surface du fluide chaud
* sfluf : surface du fluide froid
* licb  : ligne inteieure du crayon chauffant
* lecb  : ligne exteieure du crayon chauffant
* scb   : surface du crayon chauffant
* lipa lepa spa : lignes interieure, exterieure et surface du panier
* lipe lepe spe : lignes interieure, exterieure et surface du PE
* litf letf stf : lignes interieure, exterieure et surface du tdf interne
* zn       : cote axiale de la delimitation en amont de la partie construite
* zn1(n+1) : cote axiale de la delimitation en aval de la partie construite
* nmail    : nombre de maille de la tranche construite
* tabX     : tables de connectivite contenant les elements de surface
*            solides et de volume fluides.
* itab     : indice des tables
*
 
 
*** A. Objet utiles :
hnh = (zn1 - zn) / ('FLOTTANT' nmail) ;
vechau = 0. 0. hnh ;
 
*** B. Definition des elements pour une hauteur elementaire
mtfluc0 = 'VOLUME' sfluc 1 'TRAN' vechau ;
mticb0 = 'TRANSLATION' licb 1 vechau ;
mtecb0 = 'TRANSLATION' lecb 1 vechau ;
mtcb0 = 'VOLUME' scb 1 'TRAN' vechau ;
mtipa0 = 'TRANSLATION' lipa 1 vechau ;
mtepa0 = 'TRANSLATION' lepa 1 vechau ;
mtpa0 = 'VOLUME' spa 1 'TRAN' vechau ;
 
mtfluf0 = 'VOLUME' sfluf 1 'TRAN' vechau ;
mtipe0 = 'TRANSLATION' lipe 1 vechau ;
mtepe0 = 'TRANSLATION' lepe 1 vechau ;
mtpe0 = 'VOLUME' spe 1 'TRAN' vechau ;
mtitf0 = 'TRANSLATION' litf 1 vechau ;
mtetf0 = 'TRANSLATION' letf 1 vechau ;
mttf0 = 'VOLUME' stf 1 'TRAN' vechau ;
 
 
*** C. Initialisation du maillage :
   mtfluc = mtfluc0 ;
   mticb  = mticb0  ;
   mtecb  = mtecb0  ;
   mtcb   = mtcb0   ;
   mtipa  = mtipa0  ;
   mtepa  = mtepa0  ;
   mtpa   = mtpa0   ;
   mtfluf = mtfluf0 ;
   mtipe  = mtipe0  ;
   mtepe  = mtepe0  ;
   mtpe   = mtpe0   ;
   mtitf  = mtitf0  ;
   mtetf  = mtetf0  ;
   mttf   = mttf0   ;
 
*** D. Ecriture des tables de correspondances :
 
** 1ere maille :
nmtif1 = 'NBEL' mtipa ;
'REPETER' imtif1 nmtif1 ;
   lmtif1 = 'ELEM' mtecb &imtif1 ;
   lmtif3 = 'ELEM' mtipa &imtif1 ;
   lmtif5 = 'ELEM' mtepe &imtif1 ;
   lmtif7 = 'ELEM' mtitf &imtif1 ;
   tabcb . itab = lmtif1 ;
   tabpa . itab = lmtif3 ;
   tabfluc . itab = mtfluc ;
   tabpe . itab = lmtif5 ;
   tabtf . itab = lmtif7 ;
   tabfluf . itab = mtfluf ;
   itab = '+' itab 1 ;
'FIN' imtif1 ;
 
 
SI (nmail NEG 1) ;
 
'REPETER' ihaut ('-' nmail 1) ;
   vectran = 0. 0. ('*' hnh ('FLOTTANT' &ihaut)) ;
   mthe2 mtin2 mtif2 mtcb2 = PLUS mtfluc0 mticb0 mtecb0 mtcb0 vectran ;
   mtin2_p mtif2_p mtcb2_p = PLUS mtipa0 mtepa0 mtpa0 vectran         ;
   mthe4 mtin4 mtif4 mtcb4 = PLUS mtfluf0 mtipe0 mtepe0 mtpe0 vectran ;
   mtin4_p mtif4_p mtcb4_p = PLUS mtitf0 mtetf0 mttf0 vectran         ;
   nmtif2 = 'NBEL' mtif2_p ; comm changeju essai52 ;
   'REPETER' imtif2 nmtif2 ;
      lmtif2 = 'ELEM' mtif2 &imtif2 ;
      lmtif4 = 'ELEM' mtin2_p &imtif2 ;
      lmtif6 = 'ELEM' mtif4 &imtif2 ;
      lmtif8 = 'ELEM' mtin4_p &imtif2 ;
      tabcb . itab = lmtif2 ;
      tabpa . itab = lmtif4 ;
      tabfluc . itab = mthe2 ;
      tabpe . itab = lmtif6 ;
      tabtf . itab = lmtif8 ;
      tabfluf . itab = mthe4 ;
      itab = '+' itab 1 ;
   'FIN' imtif2 ;
   mtfluc = 'ET' mthe2 mtfluc ;
   mticb  = 'ET' mtin2 mticb ;
   mtecb  = 'ET' mtif2 mtecb ;
   mtcb   = 'ET' mtcb2 mtcb ;
   mtipa  = 'ET' mtin2_p mtipa ;
   mtepa  = 'ET' mtif2_p mtepa ;
   mtpa   = 'ET' mtcb2_p mtpa ;
   mtfluf = 'ET' mthe4 mtfluf ;
   mtipe  = 'ET' mtin4 mtipe ;
   mtepe  = 'ET' mtif4 mtepe ;
   mtpe   = 'ET' mtcb4 mtpe ;
   mtitf  = 'ET' mtin4_p mtitf ;
   mtetf  = 'ET' mtif4_p mtetf ;
   mttf   = 'ET' mtcb4_p mttf ;
 
'FIN' ihaut ;
 
FINSI ;
 
RESPRO mtfluc mtfluf
       mticb mtecb mtcb
       mtipa mtepa mtpa
       mtipe mtepe mtpe
       mtitf mtetf mttf
       tabcb tabpa tabpe tabtf
       tabfluc tabfluf itab ;
 
FINPRO ;
 
 
DEBPROC PLUSJU geo1*'MAILLAGE' ;
** Description
** ***********
*
* Procedure permettant de construire les maillages (lignes ou surfaces)
* a toutes les interfaces entre les parties axiales en utilisant la
* procedure 'PLUS' et la geometrie en z1:
*
 
geo2 = PLUS geo1 (0. 0. (z2 - z1)) ;
geo3 = PLUS geo1 (0. 0. (z3 - z1)) ;
geo4 = PLUS geo1 (0. 0. (z4 - z1)) ;
geo5 = PLUS geo1 (0. 0. (z5 - z1)) ;
geo6 = PLUS geo1 (0. 0. (z6 - z1)) ;
geo7 = PLUS geo1 (0. 0. (z7 - z1)) ;
geo8 = PLUS geo1 (0. 0. (z8 - z1)) ;
geo9 = PLUS geo1 (0. 0. (z9 - z1)) ;
geo10 = PLUS geo1 (0. 0. (z10 - z1)) ;
geo11 = PLUS geo1 (0. 0. (z11 -z1)) ;
geoM = PLUS geo1 (0. 0. (zM - z1)) ;
 
RESPRO geo2 geo3 geo4 geo5 geo6 geo7 geo8 geo9 geo10 geo11 geoM ;
 
FINPRO ;
 
debproc resex c*chpoint m*maillage ;
* Trace d'un sur les faces des tubes de force
* *******************************************
 
ttt = redu c m ;
mttt = extr ttt mail ;
x = prog ;
y = prog ;
z = prog ;
t = prog ;
z2 = prog ;
t2 = prog ;
n = nbno m ;
repe w n ;
 pt = mttt 'POINT' &w ;
 ptx = coor 1 pt ;
 pty = coor 2 pt ;
 ptz = coor 3 pt ;
 x = x et (prog ptx );
 y = y et (prog pty ) ;
 z = z et (prog ptz ) ;
 t9 = extr ttt scal pt ;
 t = t et (prog t9 );
 si ( et ( EGA pty 0. ) ( ptx > 0. )) ;
  z2 = z2 et ( prog ptz ) ;
  t2 = t2 et ( prog t9 ) ;
 finsi ;
fin w ;
evtnew = ORDONNER (evol manu  'Z' z2 t2) ;
finpro evtnew ;
 
 
 
 
 
DEBPROC HNAT Ttdf*'CHPOINT' Teau*'FLOTTANT' long*'FLOTTANT' ;
* Procedure permettant de calculer le coefficient d'ehange naturel
* entre le tube de force externe et la piscine
* On utilise les formules suivantes :
*   - Nu = h_nat*D/ld
*   - Nu = 0.13(Gr*Pr)**(1/3)
* Nu   : nombre de Nusselt
* long : longueur caracteristique de la convection naturelle
* ld   : conductivite de l'eau
* Gr   : nombre de Grashof
* Pr   : nombre de Prandl
* Ra   : nombre de Raileigh
* beteau : coefficient de dilatation thermique des liquides au environ
*          de 300K
* g0   : gravite* h_nat : coefficient d'ehange de convection naturelle
 
 
*Donnes du probleme :
deltaT = Ttdf - Teau ;
cnat1 = 0.13 ;
cnat2 = 1. / 3. ;
beteau = 2.E-4 ;
g0 = 9.81 ;
 
*Proprietes de l'eau :
Tfilm = KOPS (Ttdf + Teau) / 2. ;
rhoeau mueau leau cpeau = PRPEAU prt ro_o mu_o cp_o lamb_o Tfilm ;
nueau  = KOPS mueau / rhoeau ;
 
*Nombres adimensionnels :
Pr = KOPS (KOPS mueau * cpeau) / leau ;
Gr0 = g0 * beteau * ( long ** 3. ) ;
Gr = KOPS deltaT * (KOPS Gr0 / ( KOPS nueau ** 2. )) ;
Ra = KOPS Pr * Gr ;
Nu = KOPS cnat1 * ( KOPS Ra ** cnat2 ) ;
 
*Cacul du h_nat :
h_nat = KOPS Nu * (KOPS leau / long) ;
 
'RESPRO' h_nat ;
 
FINPROC ;
 
 
DEBPROC GAMMA mat*'MOT' cotx1*'FLOTTANT' $solide*'MMODEL' ;
* Procedure permettant de donner le profil axial des sources internes.
* Les parametres sont :
*    a. La composition du milieu
*    b. La position du dispositif par rapport au caisson afin de
*       deerminer l'echauffement au plan median.
*    c. Ensuite on determine le profil de l'echauffement
*
*  2. Valeur des parametres d'entree :
*    - "mat" deigne la nature du milieu. Il prend les valeurs :
*           -> "eau" pour l'eau
*           -> "alu" pour l'aluminium
*           -> "zir" pour le zirconium
*           -> "hel" pour l'heium
*           -> "azo" pour l'azote
*    - "cotx1" est la position du dispositif par rapport au caisson lors
*              de l'experience
*    - "$solide" est le modele dans lequel on calcule les sources internes dans les solides
*
*  3. Valeur des variables locales
* Eg0(x1) : sources internes au plan median a la cote x (W/g)
*           c la radioprotection qui m'a donne les mesures
* x1 : position du chariot lue sur les reglages de la manip
* z1 : position verticale en mm
* z2 : position verticale utilise par castem correspondant a z1
* Egnom   : Sources internes nominal
* Egmax   : Sources internes au plan median a z=0m
* Egcoef2 : coefficient du profil Sources internes, c un champoint
* Egch : CHPOINT donnant l'echauffemnt dans chaque maille en W/g
* Eg_SI : CHPOINT donnant l'echauffemnt dans chaque maille en W/kg
 
**Definition du facteur de Fadile :
log1 = EGA mat eau ;
log2 = EGA mat alu ;
log3 = EGA mat zir ;
log4 = EGA mat hel ;
log5 = EGA mat azo ;
 
SI (log1 OU log2 OU log4 OU log5) ;
   fad = 1. ;
FINSI ;
SI log3 ;
   fad = 1.3 ;
FINSI ;
 
Eg0 = 2.1281 * (EXP(-0.0072 * (291. - cotx1))) ;
 
 
z1 = prog -870. -700. -550. -400 -250 -175 -140 -120. -100. -80. -60.
-40. -20. 0. 20. 40. 60. 80. 100. 120. 140. 175. 250. 400. 550. 700.
860. ;
z1m = z1 * 1.e-3 ;
Egnom = prog 0. 0.065 0.13 0.25 0.675 0.84 0.861 0.895 0.925 0.948 0.97
0.985 0.995 1. 0.995 0.985 0.97 0.948 0.925 0.895 0.861 0.84 0.675 0.25
0.13 0.065 0. ;
 
 
*Changements et Interpolation :
Egcoef = Egnom * Egmax ;
 
d_soli = DOMA $solide 'MAILLAGE' ;
cotz1 = COOR 3 d_soli ; comm en m ;
 
Egcoef2 = IPOL cotz1 z1m Egcoef ; comm comparaison REFLET ;
 
**Final
Egch = KOPS (fad * Eg0) * Egcoef2 ;
Egch = NOEL $solide Egch ;
Eg_SI = Egch * 1.e3 ;
 
RESPRO Eg_SI ;
 
FINPROC ;
************************************************************************
 
 
 
 
DEBPROC FLUXNS choix*'ENTIER' lambda*'CHPOINT' T*'CHPOINT'
               $surface*'MMODEL' $solide*'MMODEL' ;
**Procedure permettant de :
**  - calculer la puissance totale traversant une interface
**  - construire un chpoint du flux (W/m2) sur cet interface
**  - choisir entre trois procedures ressemblantes
**On procede par :
**  - construction du flux avec loi de Fourrier
**  - integration sur la surface orientee par la normale sortante
 
* T est la variable temperature de la matiere du modele NS
* lambda = conductivite apparaissant dans l'equation, en argument
*         d'entree de la procedure :
*             gaine  : rvcb . inco . 'LCB'
*             panier : rvcb . inco . 'LCB_p'
* tube force interne : rvcb . inco . 'LCB_tfi'
* tube force externe : rvcb . inco . 'LCB_tfe'
*                gaz : rvcb . inco . 'LCB_gaz'
* RQ : lambda doit etre un CHPOINT scal centre !!!!!!!!!!
*
* solide  = maillage du solide
* surface = maillage de la frontiere du solide sur laquelle
*          on veut calculer le flux (Watt)
* $solide et $surface modeles s'appuyant sur ces maillages
 
**Construction des maillages a partir des modeles
d_soli = DOMA $solide 'MAILLAGE' ;
d_surf = DOMA $surface 'MAILLAGE' ;
 
SI (choix EGA 1) ;
 
*   Definition des champoints utiles a la resolution
   X Y Z = COOR d_soli ;
   R = KOPS ( KOPS X * X ) + ( KOPS Y * Y ) ;
   R = KOPS R ** 0.5 ;
   lambda = ELNO $solide lambda;
 
*   Construction du flux
   Grad_T   = KOPS T 'GRADS' $solide  ;
   Sn_Pmail = doma $surface 'XXDIAGSI' ;
*   Selon X
   Grad_Tx  = exco Grad_T Ux ;
   Flux  = KOPS Grad_Tx '*' (-1. * lambda) ;
   Fluxr = KOPS Flux * ( KOPS X / R) ;
   Fluxr = REDU Fluxr d_surf ;
*   Selon Y
   Grad_Ty  = exco Grad_T Uy ;
   Fluy  = KOPS Grad_Ty '*' (-1. * lambda) ;
   Fluyr = KOPS Fluy * ( Y / R) ;
   Fluyr = REDU Fluyr d_surf ;
 
*   Resultats :
   FluR  = KOPS Fluxr + Fluyr ;
   Wunit = KOPS FluR * Sn_Pmail ;
   Wtot  = SOMT Wunit ;
 
FINSI ;
 
 
SI (choix EGA 2) ;
 
*   Definition des champoints utiles a la resolution :
   lambda = ELNO $solide lambda; comm a evaluer en direct ;
   lambda = REDU d_surf lambda ;
   Mot1 = MOTS 'UX' 'UY' ;
   Sn_mail = DOMA $surface 'VOLUME' ;
 
*   Normale a la paroi :
   Norm0 = DOMA $solide 'NORMALEV' ;
   Norm1 = KCHT $surface VECT SOMMET Norm0 ;
 
*   Construction du flux :
   Grad_T = KOPS T 'GRADS' $solide  ;
   Grad_TS = KCHT $surface VECT SOMMET Grad_T ;
   Grad_TR = PSCA Grad_TS Norm1 Mot1 Mot1 ;
   FluR = KOPS Grad_TR '*' (-1. * lambda) ;
   FluR = NOEL $surface FluR ;
 
*   Resultats :
   Wunit = KOPS FluR * Sn_mail ;
   Wtot  = SOMT Wunit;
 
FINSI ;
 
SI (choix EGA 3) ;
 
 
*   Definition des champoints utiles a la resolution
n_surf = DOMA $surface 'MAILLAGE' ;
lambda = ELNO $solide lambda ;
lambda = REDU d_surf lambda ;
lambda = NOEL $surface lambda ;
SS1 = 'DOMA' $surface 'VOLUME' ;
 
*   Construction du flux :
   GRADC0 = 'KOPS' T 'GRAD' $solide ;
   GRADS0 = 'ELNO' $solide GRADC0 ;
   GS0    = 'KCHT' $surface 'VECT' 'SOMMET' GRADS0 ;
   GC0    = 'NOEL' $surface GS0 ;
   NORM1  = 'DOMA' $solide 'NORMALEV' ;
   NCOTE  = 'KCHT' $surface 'VECT' 'SOMMET' NORM1 ;
   NCOTE  = 'NOEL' $surface NCOTE ;
   MOT1   = 'MOTS' 'UX' 'UY' ;
 
   FLUX   = 'PSCA' GC0 NCOTE MOT1 MOT1 ;
 
   FLUX = KOPS FLUX * lambda ;
   FluR = KOPS (-1.) * FLUX ;
 
 
*   Resultats :
   Wunit = KOPS FluR * SS1 ;
   Wtot  = SOMT Wunit ;
 
FINSI ;
 
FINPROC Wtot Wunit FluR;
 
 
 
'DEBPROC'  chpobor mess1/'MOT' T*'CHPOINT' ;
*Elle permet d'afficher les valeurs minimales et maximales d'un
*champoint avec message optionnel
 
SI ( EXISTE mess1 ) ;
mess mess1 ' : ' ( MINI T ) ( MAXI T) ;
SINON ;
mess ( MINI T ) ( MAXI T) ;
FINSI ;
FINPROC ;
 
** IV.tcpu9
DEBPROC tcpu9 ;
TABTPS = TEMP 'NOEC';
tcpu = TABTPS.'TEMPS_CPU'.'INITIAL';
tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
'MESSAGE' ('CHAINE' 'tcpus=' tcpus) ;
FINPROC ;
 
'DEBP' FRAY ;
'ARGU' RVX*'TABLE' ;
 
*OPERATEUR DE RAYONNEMENT
 
RV = RVX.'EQEX' ;
* recuperation du champ de temperature
T = RV.INCO.'TCB' ;
* recuperation de la table TABR associee a l'operateur
TABR = RVX.'ARG1' ;
MAILR = TABR.'MAILLAGE' ;
* TSU : temperature de surface
TSU = 'REDU' T MAILR ;
TSU = 'EXCO' TSU 'SCAL' 'T' ;
TR = 'CHAN' 'CHAM' TSU TABR.'MODELE' 'GRAVITE' ;
* MRS : matrice de rayonnement
MRS = 'RAYN' TABR.'MODELE' TABR.'MATRICE' TR ;
MAILR = 'DOMA' TABR.'MODELENS' 'SOMMET' ;
* MRS : transformation Rigidite/Matrik
MRS = 'KOPS' 'RIMA' MRS MAILR ('MOTS' 'TCB') ;
 
AS AF = 'KOPS' 'MATRIK' ;
'FINP' AS MRS ;
 
 
 
 
 
'DEBPROC' COEFLU ma*'MAILLAGE' lch/'ENTIER' ;
* Procedure permettant de determiner les coefficients de flux lorsqu'on
* choisit d'etre a puissance lineique variable selon l'axe.
 
si ( non ( exis lch ) ) ;
cote1 = prog -700. -180.1
-180.00  -162.86  -145.71  -128.57  -111.43  -94.29  -77.14
-60.00  -42.86  -25.71   -8.57    8.57   25.71   42.86
 60.00   77.14   94.29  111.43  128.57  145.71  162.86  180.00
180.1 850. ;
 
rapP1 = prog  1.e-8 1.e-8
0.864   0.897   0.935   0.955   0.973   0.989   0.997
1.000   0.999   0.995   0.990   0.975   0.969   0.960
0.946   0.932   0.910   0.880   0.850   0.823   0.777   0.744
1.e-8 1.e-8 ;
sinon ;
 
cote1 = prog -700. -250.1
             -250.00 -226.19 -202.38 -178.57 -154.76 -130.95
             -107.14 -83.33 -59.52 -35.71 -11.90 11.90 35.71
             59.52 83.33 107.14 130.95 154.76 178.57 202.38
             226.19 250.00 250.1 850. ;
rapP1 = prog 1.e-8 1.e-8
             0.722 0.747 0.812 0.867 0.914 0.952 0.977
             0.996 1.000 0.997 0.992 0.972 0.966 0.947 0.924
             0.888 0.847 0.801 0.748 0.689 0.622 0.548
             1.e-8 1.e-8 ;
 
 
finsi ;
*Changement d'unite : on passe en metre
cote3 = cote1 * 1e-3 ;
 
 
** yyv1 est un CHAMPOINT representant le profil de flux
xv1 = ( COOR 3 ma ) ;
yyv1 =  IPOL xv1 cote3 rapP1 ;
 
'RESPRO' yyv1 ;
'FINPROC' ;
 
 
DEBPROC EBULL ;
*   Procedure permettant de determiner si il y a ou non ebullition locale
* sur la gaine du crayon chauffant. Si oui, cette procedure donne
* numeriquement et graphiquement la zone exacte ou apparait l'ebullition
* nuclee.
*   Les reultats obtenues dependent fortement du maillage, en particulier
* du nombre de maille selon z : nhaut.
*   On procede de la maniere suivante
*     - Determination de la temperature de paroi a saturation :
* a l'aide de la formule de Jens et Lottes
*     - Determination des coefficiant indiquant la zone d'ebullition :
* c_mono : coefficient d'echange monophasique =| 0 si ebullition nuclee
*                                              | 1 sinon
* c_di   : coefficient d'echange diphasique =| 1 si ebullition nuclee
*                                            | 0 sinon
*     - Determination de la cote d'ebullition
* ebull1 =| VRAI si il y a ebullition
*         | FAUX si on reste en monophasique
*
* L'indice "_N" deigne le champoint correspondant aux noeuds du
* maillage
 
***Remise a niveau de f0 ;
cf0 = f0 MASQ 'SUPERIEUR' 0. ;
f0 = KOPS f0 * cf0 ;
***Determination de la temperature de paroi a saturation :
Volif = DOMA $mtif 'VOLUME' ;
 phi0 = KOPS f0 / Volif;
 phi0 = 1.e-4 * phi0 ;
 
expo1 = EXP ( -1. * PsatB / 62.) ;
DTsat = 7.91 * expo1 *(phi0 ** 0.25) ;
Tpsat = Tsat + DTsat ;
DTsat_N = ELNO $mtif DTsat ; comm modif essai14 ;
Tpsat_N = Tsat + DTsat_N ; comm modif essai14 ;
 
***Determination des coeff indiquant la zone d'ebullition
Tparoi = NOEL $mtif (rvi . 'TCB') ;
c_mono = Tparoi MASQ 'INFERIEUR' Tpsat ;
c_diph = Tparoi MASQ 'SUPERIEUR' Tpsat ;
 c_app = Tparoi MASQ 'SUPERIEUR' 'SOMME' Tpsat ;
 
Tparoi_N = REDU (rvi . 'TCB') mtif ; comm modif essai14 ;
c_mono_N = Tparoi_N MASQ 'INFERIEUR' Tpsat_N ; comm modif essai14 ;
c_diph_N = Tparoi_N MASQ 'SUPERIEUR' Tpsat_N ; comm modif essai14 ;
 
SI (c_app EGA 0) ;
   ebulli = FAUX ;
   RESPRO ebulli c_mono c_diph c_mono_N c_diph_N ;
 
SINON ;
   ebulli = VRAI ;
***Determination de la cote d'ebullition
 
*      1. Numeriquement :
   Z = COOR 3 mtif ;
   Zeb = Z * c_diph_N ; comm modif essai14 ;
   mtif_pt = CHAN POI1 mtif ;
   nbno1 = nbno mtif ;
   listZeb = prog ;
 
   'REPETER' boucju nbno1 ;
      pt1 = mtif_pt point &boucju  ;
      valZeb = extr Zeb 'SCAL' pt1    ;
      SI (valZeb > 0) ;
         listZeb = listZeb ET (prog valZeb) ;
      FINSI ;
   FIN boucju ;
 
   ZebMin = MINI listZeb ;
   ZebMax = MAXI listZeb ;
   MESS 'ebullition nuclee apparait entre :' ZebMin 'et' ZebMax 'm' ;
 
***Hauteur des mailles dans les differentes parties :
   hautmA = (z2 - z1) / nhautA ;
   hautmB = (z3 - z2) / nhautB ;
   hautmC = (z4 - z3) / nhautC ;
   hautmD = (z5 - z4) / nhautD ;
   hautmE = (zM - z5) / nhautE ;
 
   Mess 'taille d une maille partie A' hautmA ;
   Mess 'taille d une maille partie B' hautmB ;
   Mess 'taille d une maille partie C' hautmC ;
   Mess 'taille d une maille partie D' hautmD ;
   Mess 'taille d une maille partie E' hautmE ;
 
*      2. Graphiquement
*   Evc_diph = EVOL ROUG CHPO c_diph_N ligif; comm modif essai14 ;
*   DESS Evc_diph 'TITR' 'Zone ebullition nuclee' 'TITX' 'Hauteur' ;
 
   RESPRO ebulli c_mono c_diph c_mono_N c_diph_N ;
FINSI ;
FINPROC ;
 
************************************************************************
**************************FIN DES PROCEDURES****************************
************************************************************************
*opti debu 1;
 
 
 
discr    = 'LINE' ;
 
***graffic  = vrai ; comm trace fichier ;
graffic  = faux ; comm trace ecran ;
 
interact = VRAI   ;
interact = faux   ; comm execution sans arret ;
 
rayo     = VRAI   ;
 
egasol = VRAI ;
egaliq = VRAI ;
 
verbose  = VRAI   ;
precis   = VRAI   ;
 
complet = vrai ;
complet = faux ;
 
si ( non complet ) ;
 rayo = faux ;
finsi ;
 
 
si interact  ;
 graph = vrai ;
sinon ;
 graph = faux ;
finsi ;
* graph = vrai ; opti 'TRAC' 'PSC' 'EPTR' 6 'POTR' 'HELVETICA_16';
 
graph0    = graph;
graffic0  = graffic ;
interac0 = interact ;
rayo0     = rayo ;
egasol0 = egasol ;
egaliq0 = egaliq ;
 
verbose0  = verbose ;
precis0   = precis ;
comple0   = complet ;
 
 
** II.Parametres geometriques et numeriques
*      II.1.Description
*
* nb_bouc : nombre maximal d'iterations de la boucle
* ne      : nombre de maille selon r
* np      : nombre de maiile selon theta
*
* prhe   : pression de l'helium dans la cavite du tube de force, en bars
* qmhe   : debit massique des fluides chaud et froid, en kg/s
* theec  : |temperature, en Celsius,|d'entree du fluide chaud
*                                  |de sortie du fluide froid
*          |imposee en realite par les systemes de regulation d'ISA1
* tpisci : temperature de la piscine, en Celsius
* plin1  : |puissance lineique delivree par le crayon chauffant
*          |en W.cm-1
* PsatB  : pression constante des fluides et donc de saturation, en Bars
* CK      : permet de passer des degre Celsius en Kelvin.
* posix   : position du dispositif par rapport au caisson
*         : valeur que l'on lit, sur les ecrans de controle de la manip
*           ISA1, en mm
 
 
*      II.2.Valeurs
nb_bouc = 40 ;
SI rayo ;
      ne = 3 ;
      np = 22 ; comm suffisamment grands sinon rayo marche pas ;
      nral = 10 ; comm discretisation radiale des rallonges ;
SINON ;
      ne = 3 ;
      np = 10;
      nral = 10 ; comm discretisation radiale des rallonges ;
FINSI ;
 
* Loi de convection forcee
A_colb = 0.023 ;
B_colb = 0.8D0 ;
C_colb = ('/' 1.D0 3.D0) ;
 
prhe   = 5.D0   ;
 
**Donnees variables du calcul
**--------------------------
 
qmhe   = 0.1956 ;
theec  = 252.   ;
tpisci = 40.    ;
plin1  = 480.   ;
PsatB  = 146.4  ;
Posix  = 198.   ;
Egmax  = 1.089 ;
***Coefficient echange externe piscine
***valeur positive�: coeff constant sur la hauteur
***valeur negative�: convection naturelle et calcul auto
Hcnczt = -1000. ;
diamex0 = 9.5D-3 ;
epai0 = 1.D-3 ;
 
lat00 = 1045713.55 ; comm pression au dela des tables Cast3M ;
 
mess 'valeurs lues : ' ;
mess num0 qmhe theec plin1 Egmax
   posix psatb diamex0 epai0 hcnczt ;
 
si ( hcnczt &lt;EG 0 ) ;
 varcna = vrai ;
sinon ;
 varcna = faux ;
 hcncst = hcnczt ;
finsi ;
 
SI varcna ;
 mess 'Convection naturelle avec la piscine' ;
SINON ;
 mess 'Echange constant avec la piscine : ' hcncst ' W/m2.K' ;
FINSI ;
 
varcna0   = varcna ;
hcncst0   = hcncst ;
 
ppa   = PsatB * 1.e5 ; comm pression en Pascal ;
CK     = 273.15 ;
 
 
prhe0   = prhe ;
qmhe0   = qmhe ;
theec0  = theec ;
theeK0  = theeK ;
tpisci0 = tpisci ;
plin10  = plin1 ;
psatb0  = psatb ;
ppa0    = ppa ;
posix0 = posix ;
Egmax0 = Egmax ;
 
choiju = 2 ; comm choix de la procedure calculant les flux ;
 
* PARAMETRAGE
************************************************************************
'TEMPS' 'ZERO' ;
'OPTION' 'DIME' 3 ;
'OPTION' 'ISOV' 'SULI' ;
nbisov = 15 ;
si graffic ;
  'OPTION' 'TRAC' 'PSC' ;
*  'OPTION' 'ECHO' 0 ;
'SINON' ;
  'OPTION' 'TRAC' 'X' ;
'FINSI' ;
 
titglob = 'CHAINE' ' Isabelle1' ;
 
**ESSAI30:**
'SI' ('EGA' discr 'LINE') ;
   'OPTION' 'DIME' 3 'ELEM' 'CUB8' ;
'SINON' ;
   'OPTION' 'DIME' 3 'ELEM' 'CU20' ;
'FINSI' ;
**
 
************************************************************************
 
 
* MAILLAGE
************************************************************************
** I.Geometrie
*      I.1.Parametres initiaux
*
*             I.1.1.Description
*
* partie A => rallonge infeieure,        entre z1 et z2 ;
* partie B => bouchon infeieur  ,        entre z2 et z3 ;
* partie C => crayon             ,        entre z3 et z4 ;
* partie D => bouchon supeieur  ,        entre z4 et z5 ;
* partie E => rallonge superieure,        entre z5 et z6 ;
* partie F => changement section tdf ext, entre z6 et z7 ;
* partie G => changement section tdf int, entre z7 et z8 ;
* partie H => jusqu'au fin du panier,     entre z8 et z9 ;
* partie I => depuis la fin du panier,    entre z9 et z10 ;
* partie J => changement section fluid ch,entre z10 et z11 ;
* partie K => plus de rallonge sup,       entre z11 zt zM ;
*
*
*             I.1.1.Valeurs numeriques
** COTES INVARIANTES
* z1 : coordonne axiale du bas du dispositif modelise,  (m).
* zM : coordonne axiale du haut du dispositif modelise, (m).
* o : point origine
* o1 : point de l'axe
 
z1 = -695D-3  ;
z6 = 539.D-3  ;
z7 = 569.D-3 ;
z8 = 597.D-3 ;
z9 = 630.D-3 ;;
zM = 850D-3   ;
o = 0. 0. 0.  ;
o1 = 0. 0. 1. ;
 
** COTES VARIANTES
z2 = -197.4D-3 ;
*** Cotes du crayon
z3 = -180.D-3 ;
z4 = 180.D-3 ;
z5 = 248.1D-3 ;
z10 = 743.D-3 ;
z11 = z10 + 27.D-3 ;
 
** Nombre de maille selon la hauteur
nhautB = 1 ;
nhautD = 1 ;
nhautF = 1 ;
nhautG = 1 ;
nhautH = 1 ;
nhautJ = 1 ;
nhautK = 1 ;
si complet ;
 nhautA = 6 ;
 nhautE = 4 ;
 nhautC = 20 ;
 nhautI = 2 ;
sinon ;
 nhautA = 3 ;
 nhautE = 2 ;
 nhautC = 5 ;
 nhautI = 1 ;
finsi ;
 
 
 
 
*** LE CRAYON :
* partie A => rallonge inferieure (Zircaloy4)
* partie B => bouchon inferieur (vide)
* partie C => longueur chauffante du crayon
* partie D => bouchon superieur (vide)
* parties E,F,G,H,I => rallonge superieure (Zircaloy4)
* diamext : diametre exterieur de la gaine entourant le crayon chauffant
* epai    : epaisseur de la gaine entourant le crayon chauffant
diamext = diamex0 ;
epai = epai0 ;
****diamext = 9.5D-3  ;
****epai    = 1.D-3   ;
 
*** LE PANIER :
* pani : diametre interieur du panier (m)
* pane : diametre exterieur du panier (m)
pani = 19.6D-3 ;
pane = 21.0D-3 ;
 
*** LE FLUIDE CHAUD ASCENDANT:
* spass : surface de passage du fluide chaud, i-e l'ecoulement interne
* cotp  : cote du carre representant la base du maillage du fluide chaud
* dhy   : diametre hydraulique du fluide chaud
* qhenom   : vitesse massique du fluide chaud (kg.m-2.s-1)
 
spass   = (pi * ((pani / 2.) ** 2.)) - (pi * ((diamext / 2.) ** 2.)) ;
cotp    = '**' spass 0.5D0 ;
dhy     = 4 * spass / (( pi * diamext ) + ( pi * pani)) ;
qhenom  = qmhe / spass ;
 
spassK  = (pi * ((pani / 2.) ** 2.)) ;
cotpK   = '**' spassK 0.5D0 ;
dhyK    = 4 * spassK / ( pi * pani) ;
qhenomK = qmhe / spassK ;
 
 
*** LE PORTE ECHANTILLON (PE) :
* pei : diametre interieur du PE (m)
* pee : diametre exterieur du PE (m)
*****pei = 23.0D-3 ;
pei = 22.0D-3 ; comm DRSN/SIREN/LECSI/DO 16 du 11/02/08 ;
pee = 25.0D-3 ;
 
*** LE TUBE DE FORCE INTERNE :
* partie A,B,C,D,E => 1ere partie tdf interne (Zircaloy4)
* partie F => changement de section du tdf externe (Zircaloy4)
* partie G => changement de section du tdf interne (Zircaloy4)
* partie H,I,J,K => 2eme partie tdf interne et externe (Zircaloy4)
* tfintiA : diametre interne du tube de force interne partie A -> E
* tfinteA : diametre externe du tube de force interne partie A -> E
* tfintiH : diametre interne du tube de force interne partie H -> J
* tfinteH : diametre externe du tube de force interne partie H -> J
 
tfintiA  = 46.D-3   ;
tfinteA  = 59.0D-3  ;
tfintiH  = 53.05D-3 ;
tfinteH  = 66.50D-3 ;
 
*** LE FLUIDE FROID DESCENDANT:
* spassN  : surface de passage dans la partie N
* cotpN   : cote du carre representant la base du maillage partie N
* dhyN    : diametre hydraulique de la partie N
* qhenomN : vitesse massique de la partie N (kg.m-2.s-1)
 
spassA  = (pi * ((tfintiA / 2.) ** 2.)) - (pi * ((pee / 2.) ** 2.)) ;
cotpA   = '**' spassA 0.5D0 ;
dhyA    = 4 * spassA / (( pi * pee ) + ( pi * tfintiA)) ;
qhenomA = qmhe / spassA ;
 
spassH  = (pi * ((tfintiH / 2.) ** 2.)) - (pi * ((pee / 2.) ** 2.)) ;
cotpH   = '**' spassH 0.5D0 ;
dhyH    = 4 * spassH / (( pi * pee ) + ( pi * tfintiH)) ;
qhenomH = qmhe / spassH ;
 
 
*** LE TUBE DE FORCE EXTERNE :
* parties A,B,C,D,E => 1ere partie tdf externe (Zircaloy4)
* partie F => changement de section du tdf externe (Zircaloy4)
* partie G => changement de section du tdf interne (Zircaloy4)
* partie H,I,J,K => 2eme partie tdf externe (Zircaloy4)
* tfextiA : diametre interne du tube de force externe partie A -> E
* tfexteA : diametre externe du tube de force externe partie A -> E
* tfextiG : diametre interne du tube de force externe partie G -> J
* tfexteG : diametre externe du tube de force externe partie G -> J
 
*****tfextiA  = 59.9D-3 ;
*****tfexteA  = 70.1D-3 ;
tfextiA  = 59.8D-3 ; comm DRSN/SIREN/LECSI/DO 16 du 11/02/08 ;
tfexteA  = 70.D-3 ; comm DRSN/SIREN/LECSI/DO 16 du 11/02/08 ;
tfextiG  = 67.5D-3 ;
tfexteG  = 78.5D-3 ;
 
 
*      I.2.Parametres calcules
*             I.2.1.Description
* theek    : conversion en Kelvin de theec
* tpiscik  : conversion en Kelvin de tpisci
* plin1    : conversion en W/m.
* Psat : conversion en Pascals de PsatB
* Tsat : temperature de saturation a la pression Psat
 
*             I.2.2.Valeurs
 
theek   = theec + CK ;
tpiscik = tpisci + CK ;
plin1   = '*' plin1 100.D0 ;
Psat = PsatB * 1.e5 ;
Tsat = CK + ( h2o_tsat Psat );
** IV.Points, lignes, surfaces
*      IV.1. Crayon Chauffant
 
* A. GAINE :
rext = diamext / 2. ;
rint = rext - epai ;
mrext = -1. * rext ;
mrint = -1. * rint ;
 
**Points de l'interface 1 :
p11e = rext 0. z1 ;
p12e = 0. rext z1 ;
p13e = mrext 0. z1 ;
p14e = 0. mrext z1 ;
p11i = rint 0. z1 ;
p12i = 0. rint z1 ;
p13i = mrint 0. z1 ;
p14i = 0. mrint z1 ;
 
*** Lignes de toutes les interfaces :
l11e = 'CER3' np p11e p12e p13e ;
l12e = 'CER3' np p13e p14e p11e ;
lt1e = 'ET' l11e l12e ;
l11i = 'CER3' np p11i p12i p13i ;
l12i = 'CER3' np p13i p14i p11i ;
lt1i = 'ET' l11i l12i ;
 
lt2e lt3e lt4e lt5e lt6e lt7e lt8e lt9e lt10e lt11e noMe = PLUSJU lt1e ;
lt2i lt3i lt4i lt5i lt6i lt7i lt8i lt9i lt10i lt11i noMi = PLUSJU lt1i ;
 
*** Bases des volumes :
stcb1 = 'REGLER' lt1i ne lt1e ;
stcb2 stcb3 stcb4 stcb5 stcb6 stcb7 stcb8 stcb9 stcb10 stcb11 noM
 = PLUSJU stcb1 ; comm meme si stcb11 n existe pas ;
 
ELIM 1.e-8 (lt2e ET lt3e ET lt4e ET lt5e ET lt6e ET lt7e ET lt8e ET
            lt9e ET lt10e ET stcb1 ET stcb2 ET stcb3 ET stcb4
            ET stcb5 ET stcb6 ET stcb7 ET stcb8 ET stcb9 ET stcb10) ;
 
* B. LES RALLONGES :
rintr = rint / 10. ;
mrintr = -1. * rintr ;
 
**Points :
p11r = rintr 0. z1 ;
p12r = 0. rintr z1 ;
p13r = mrintr 0. z1 ;
p14r = 0. mrintr z1 ;
 
*** Lignes et Surfaces : interfaces 1->10 :
l11r = 'CER3' np p11r p12r p13r ;
l12r = 'CER3' np p13r p14r p11r ;
lt1r = 'ET' l11r l12r ;
 
stcb11r = 'REGLER' lt1r nral lt1i ;
stcb12r = SURF lt1r 'PLANE' ;
stcb1r =  stcb11r et stcb12r ;
stcb2r stcb3r stcb4r stcb5r stcb6r stcb7r stcb8r stcb9r stcb10r no11 noM
 = PLUSJU stcb1r ;
 
 
*      IV.2. Fluide chaud     ESSAI30 :
 
* 1. Section 1 : A->I
*RQ : La geometrie est un carre de surface egale a spass
*     Au lieu de centrer le canal du fluide chaud en (0;0) on va le
*     centrer en (0.07;0).
 
*** Points de l'interface 1 :
cotp2 = '/' cotp 2.D0 ;
mcotp2 = '*' cotp2 -1.D0 ;
p11 = (70.e-3 + cotp2) cotp2 z1 ;
p12 = (70.e-3 + mcotp2) cotp2 z1 ;
p13 = (70.e-3 + mcotp2) mcotp2 z1 ;
p14 = (70.e-3 +cotp2)  mcotp2 z1 ;
 
*** Surfaces :
sthe1 = 'MANUEL' 'QUA4' p11 p12 p13 p14;
sthe2 sthe3 sthe4 sthe5 sthe6 sthe7 sthe8 sthe9 sthe10 no11 noM
 = PLUSJU sthe1 ;
 
* 2. Section 2 :
*La geometrie est un carre de surface egale a spassK
cotp2K = '/' cotpK 2.D0 ;
mcotp2K = '*' cotp2K -1.D0 ;
p111 = (70.e-3 + cotp2K) cotp2K z11 ;
p112 = (70.e-3 + mcotp2K) cotp2K z11 ;
p113 = (70.e-3 + mcotp2K) mcotp2K z11 ;
p114 = (70.e-3 +cotp2K)  mcotp2K z11 ;
sthe11 = 'MANUEL' 'QUA4' p111 p112 p113 p114 ;
stheM = PLUS sthe11 (0. 0. (zM - z11)) ;
 
'SI' ('NON' ('EGA' discr 'LINE')) ;
   sthe1 = 'CHANGER' 'QUADRATIQUE' sthe1 ;
   sthe2 = 'CHANGER' 'QUADRATIQUE' sthe2 ;
   sthe3 = 'CHANGER' 'QUADRATIQUE' sthe3 ;
   sthe4 = 'CHANGER' 'QUADRATIQUE' sthe4 ;
   sthe5 = 'CHANGER' 'QUADRATIQUE' sthe5 ;
   sthe6 = 'CHANGER' 'QUADRATIQUE' sthe6 ;
   sthe7 = 'CHANGER' 'QUADRATIQUE' sthe7 ;
   sthe8 = 'CHANGER' 'QUADRATIQUE' sthe8 ;
   sthe9 = 'CHANGER' 'QUADRATIQUE' sthe9 ;
   sthe10 = 'CHANGER' 'QUADRATIQUE' sthe10 ;
   sthe11 = 'CHANGER' 'QUADRATIQUE' sthe11 ;
   stheM = 'CHANGER' 'QUADRATIQUE' stheM ;
'FINSI' ;
 
*      IV.3. Panier(interne)     ESSAI30 :
 
rint1 = ( pani / 2.D0 ) ;
rext1 = ( pane / 2.D0 ) ;
mrext1 = '*' rext1 -1.D0 ;
mrint1 = '*' rint1 -1.D0 ;
 
**Points de l'interface 1 :
q11e = rext1 0. z1    ;
q12e = 0. rext1 z1  ;
q13e = mrext1 0. z1 ;
q14e = 0. mrext1 z1 ;
q11i = rint1 0. z1  ;
q12i = 0. rint1 z1  ;
q13i = mrint1 0. z1 ;
q14i = 0. mrint1 z1 ;
 
*** Lignes de toutes les interfaces :
l11e_p = 'CER3' np q11e q12e q13e ;
l12e_p = 'CER3' np q13e q14e q11e ;
lt1e_p = 'ET' l11e_p l12e_p ;
l11i_p = 'CER3' np q11i q12i q13i ;
l12i_p = 'CER3' np q13i q14i q11i ;
lt1i_p = 'ET' l11i_p l12i_p ;
 
lt2e_p lt3e_p lt4e_p lt5e_p lt6e_p lt7e_p lt8e_p lt9e_p lt10e_p lt11e_p
 ltMe_p = PLUSJU lt1e_p ;
lt2i_p lt3i_p lt4i_p lt5i_p lt6i_p lt7i_p lt8i_p lt9i_p lt10i_p lt11i_p
 ltMi_p = PLUSJU lt1i_p ;
 
*** Surfaces de toutes les interfaces :
stcb_p1 = 'REGLER' lt1i_p ne lt1e_p ;
stcb_p2 stcb_p3 stcb_p4 stcb_p5 stcb_p6 stcb_p7 stcb_p8 stcb_p9
 stcb_p10 stcb_p11 stcb_pM = PLUSJU stcb_p1 ;
 
ELIM 1.e-8 (lt2e_p ET lt3e_p ET lt4e_p ET lt5e_p ET lt6e_p ET lt7e_p ET
            lt8e_p ET lt9e_p ET lt10e_p ET lt11e_p ET ltMe_p ET lt2i_p
            ET lt3i_p ET lt4i_p ET lt5i_p ET lt6i_p ET lt7i_p ET lt8i_p
            ET lt9i_p ET lt10i_p ET lt11i_p ET ltMi_p ET stcb_p1
            ET stcb_p2 ET stcb_p3 ET stcb_p4 ET stcb_p5 ET stcb_p6 ET
            stcb_p7 ET stcb_p8 ET stcb_p9 ET stcb_p10 ET stcb_p11 ET
            stcb_pM) ;
 
 
*      IV.4. Porte echantillon     ESSAI32 :
 
rpeint = (pei / 2.D0 ) ;
rpeext = ( pee / 2.D0 ) ;
mrpeint = '*' rpeint -1.D0 ;
mrpeext = '*' rpeext -1.D0 ;
 
**Points de l'interface 1 :
qpe11e = rpeext 0. z1  ;
qpe12e = 0. rpeext z1  ;
qpe13e = mrpeext 0. z1 ;
qpe14e = 0. mrpeext z1 ;
qpe11i = rpeint 0. z1  ;
qpe12i = 0. rpeint z1  ;
qpe13i = mrpeint 0. z1 ;
qpe14i = 0. mrpeint z1 ;
 
*** Lignes de toutes les interfaces :
lpe11e = 'CER3' np qpe11e qpe12e qpe13e ;
lpe12e = 'CER3' np qpe13e qpe14e qpe11e ;
lpet1e = 'ET' lpe11e lpe12e ;
lpe11i = 'CER3' np qpe11i qpe12i qpe13i ;
lpe12i = 'CER3' np qpe13i qpe14i qpe11i ;
lpet1i = 'ET' lpe11i lpe12i ;
 
lpet2e lpet3e lpet4e lpet5e lpet6e lpet7e lpet8e lpet9e lpet10e lpet11e
 lpetMe = PLUSJU lpet1e ;
lpet2i lpet3i lpet4i lpet5i lpet6i lpet7i lpet8i lpet9i lpet10i lpet11i
 lpetMi = PLUSJU lpet1i ;
 
*** Surfaces de toutes les interfaces :
stpe1 = 'REGLER' lpet1i ne lpet1e ;
stpe2 stpe3 stpe4 stpe5 stpe6 stpe7 stpe8 stpe9 stpe10 stpe11 stpeM
 = PLUSJU stpe1 ;
 
ELIM 1.e-8 (lpet2e ET lpet3e ET lpet4e ET lpet5e ET lpet6e ET lpet7e ET
            lpet8e ET lpet9e ET lpet10e ET lpet11e ET lpetMe ET stpe1
            ET stpe2 ET stpe3 ET stpe4 ET stpe5 ET stpe6 ET stpe7 ET
            stpe8 ET stpe9 ET stpe10 ET stpe11 ET stpeM) ;
 
 
*      IV.5. Lame d'eau stagnante : ESSAI51 :
 
st_lam1 = 'REGLER' lt1e_p ne lpet1i ;
st_lam2 st_lam3 st_lam4 st_lam5 st_lam6 st_lam7 st_lam8 st_lam9 st_lam10
 st_lam11 st_lamM = PLUSJU st_lam1 ;
 
 
*      IV.6. Fluide froid     ESSAI32 :
 
* 1. Section 1 :
*La geometrie est un carre de surface egale a spassA
*     Au lieu de centrer le canal du fluide chaud en (0;0) on va le
*     centrer en (0.;0.08).
 
cotp12 = cotpA / 2. ;
mcotp12 = -1.D0 * cotp12 ;
q11 = cotp12 (80.e-3 + cotp12) z1 ;
q12 = mcotp12 (80.e-3 + cotp12) z1 ;
q13 = mcotp12 (80.e-3 + mcotp12) z1 ;
q14 = cotp12  (80.e-3 + mcotp12) z1 ;
 
sthe_p1 = 'MANUEL' 'QUA4' q11 q12 q13 q14 ;
sthe_p2 sthe_p3 sthe_p4 sthe_p5 sthe_p6 sthe_p7 no8 no9 no10 no11 noM
 = PLUSJU sthe_p1 ;
 
 
* 2. Section 2 :
*La geometrie est un carre de surface egale a spassH
cotp22 = cotpH / 2. ;
mcotp22 = -1.D0 * cotp22 ;
q81 = cotp22 (80.e-3 + cotp22)  z8 ;
q82 = mcotp22 (80.e-3 + cotp22) z8 ;
q83 = mcotp22 (80.e-3 + mcotp22) z8 ;
q84 = cotp22 (80.e-3 + mcotp22) z8 ;
sthe_p8 = 'MANUEL' 'QUA4' q81 q82 q83 q84 ;
sthe_p9 = PLUS sthe_p8 (0. 0. (z9 - z8)) ;
sthe_p10 = PLUS sthe_p8 (0. 0. (z10 - z8)) ;
sthe_p11 = PLUS sthe_p8 (0. 0. (z11 - z8)) ;
sthe_pM = PLUS sthe_p8 (0. 0. (zM - z8)) ;
 
'SI' ('NON' ('EGA' discr 'LINE')) ;
   sthe_p1 = 'CHANGER' 'QUADRATIQUE' sthe_p1 ;
   sthe_p2 = 'CHANGER' 'QUADRATIQUE' sthe_p2 ;
   sthe_p3 = 'CHANGER' 'QUADRATIQUE' sthe_p3 ;
   sthe_p4 = 'CHANGER' 'QUADRATIQUE' sthe_p4 ;
   sthe_p5 = 'CHANGER' 'QUADRATIQUE' sthe_p5 ;
   sthe_p6 = 'CHANGER' 'QUADRATIQUE' sthe_p6 ;
   sthe_p7 = 'CHANGER' 'QUADRATIQUE' sthe_p7 ;
   sthe_p8 = 'CHANGER' 'QUADRATIQUE' sthe_p8 ;
   sthe_p9 = 'CHANGER' 'QUADRATIQUE' sthe_p9 ;
   sthe_p10 = 'CHANGER' 'QUADRATIQUE' sthe_p10 ;
   sthe_p11 = 'CHANGER' 'QUADRATIQUE' sthe_p11 ;
   sthe_pM = 'CHANGER' 'QUADRATIQUE' sthe_pM ;
'FINSI' ;
 
 
*      IV.7. Tube de force interne      ESSAI34 :
 
rtiinA = (tfintiA / 2.D0) ;
rtiexA = (tfinteA / 2.D0) ;
mrtiexA = '*' rtiexA -1.D0 ;
mrtiinA = '*' rtiinA -1.D0 ;
rtiinH = (tfintiH / 2.D0) ;
rtiexH = (tfinteH / 2.D0) ;
mrtiexH = '*' rtiexH -1.D0 ;
mrtiinH = '*' rtiinH -1.D0 ;
 
**Points de l'interface 1 :
qti11e = rtiexA 0. z1  ;
qti12e = 0. rtiexA z1  ;
qti13e = mrtiexA 0. z1 ;
qti14e = 0. mrtiexA z1 ;
qti11i = rtiinA 0. z1  ;
qti12i = 0. rtiinA z1  ;
qti13i = mrtiinA 0. z1 ;
qti14i = 0. mrtiinA z1 ;
 
** Points de l'interface 7 :
qti71e = rtiexA 0. z7  ;
qti71i = rtiinA 0. z7  ;
 
**Points de l'interface 8 :
qti81e = rtiexH 0. z8  ;
qti82e = 0. rtiexH z8  ;
qti83e = mrtiexH 0. z8 ;
qti84e = 0. mrtiexH z8 ;
qti81i = rtiinH 0. z8  ;
qti82i = 0. rtiinH z8  ;
qti83i = mrtiinH 0. z8 ;
qti84i = 0. mrtiinH z8 ;
 
* Points supplementaires pour la cavite de gaz :
qti21e = rtiexA 0. z2  ;
qti31e = rtiexA 0. z3  ;
qti41e = rtiexA 0. z4  ;
qti51e = rtiexA 0. z5  ;
qti61e = rtiexA 0. z6  ;
qti91e = rtiexH 0. z9  ;
qti101e = rtiexH 0. z10  ;
qti111e = rtiexH 0. z11  ;
qtiM1e = rtiexH 0. zM  ;
 
*** Lignes des interfaces 1 -> 7 :
lti11e_p = 'CER3' np qti11e qti12e qti13e ;
lti12e_p = 'CER3' np qti13e qti14e qti11e ;
ltit1e_p = 'ET' lti11e_p lti12e_p ;
lti11i_p = 'CER3' np qti11i qti12i qti13i ;
lti12i_p = 'CER3' np qti13i qti14i qti11i ;
ltit1i_p = 'ET' lti11i_p lti12i_p ;
 
ltit2e_p ltit3e_p ltit4e_p ltit5e_p ltit6e_p ltit7e_p no8 no9 no10 no11
 noM = PLUSJU ltit1e_p ;
ltit2i_p ltit3i_p ltit4i_p ltit5i_p ltit6i_p ltit7i_p no8 no9 no10 no11
 noM = PLUSJU ltit1i_p ;
 
*** Lignes des interfaces 8 -> 11 :
lti81e_p = 'CER3' np qti81e qti82e qti83e ;
lti82e_p = 'CER3' np qti83e qti84e qti81e ;
ltit8e_p = 'ET' lti81e_p lti82e_p ;
lti81i_p = 'CER3' np qti81i qti82i qti83i ;
lti82i_p = 'CER3' np qti83i qti84i qti81i ;
ltit8i_p = 'ET' lti81i_p lti82i_p ;
 
ltit9e_p  = PLUS ltit8e_p (0. 0. (z9 - z8))  ;
ltit10e_ = PLUS ltit8e_p (0. 0. (z10 - z8)) ;
ltit11e_ = PLUS ltit8e_p (0. 0. (z11 - z8)) ;
ltit9i_p = PLUS ltit8i_p (0. 0. (z9 - z8)) ;
ltit10i_ = PLUS ltit8i_p (0. 0. (z10 - z8)) ;
ltit11i_ = PLUS ltit8i_p (0. 0. (z11 - z8)) ;
 
*** Surfaces de toutes les interfaces :
stfint1 = 'REGLER' ltit1i_p ne ltit1e_p ;
stfint2 = 'REGLER' ltit2i_p ne ltit2e_p ;
stfint3 = 'REGLER' ltit3i_p ne ltit3e_p ;
stfint4 = 'REGLER' ltit4i_p ne ltit4e_p ;
stfint5 = 'REGLER' ltit5i_p ne ltit5e_p ;
stfint6 = 'REGLER' ltit6i_p ne ltit6e_p ;
stfint7 = 'REGLER' ltit7i_p ne ltit7e_p ;
stfint8 = 'REGLER' ltit8i_p ne ltit8e_p ;
stfint9 = 'REGLER' ltit9i_p ne ltit9e_p ;
stfint10 = 'REGLER' ltit10i_ ne ltit10e_ ;
stfint11 = 'REGLER' ltit11i_ ne ltit11e_ ;
stfintM = PLUS stfint11 (0. 0. (zM - z11)) ;
 
elim 1.e-8 (stfint2 ET stfint3 ET stfint4 ET stfint5 ET stfint6 ET
            stfint7 ET stfint9 ET stfint10 ET stfint11 ET stfintM ET
            qti21e ET qti31e ET qti41e ET qti51e ET qti61e ET qti71e ET
            qti91e ET qti101e ET qti111e ET qtiM1e) ;
 
 
*      IV.8. Tube de force externe
 
rteinA = ( tfextiA / 2.D0 ) ;
rteexA = ( tfexteA / 2.D0 ) ;
mrteexA = '*' rteexA -1.D0 ;
mrteinA = '*' rteinA -1.D0 ;
rteinG = ( tfextiG / 2.D0 ) ;
rteexG = ( tfexteG / 2.D0 ) ;
mrteexG = '*' rteexG -1.D0 ;
mrteinG = '*' rteinG -1.D0 ;
 
**Points de l'interface 1 :
qte11e = rteexA 0. z1 ;
qte12e = 0. rteexA z1 ;
qte13e = mrteexA 0. z1 ;
qte14e = 0. mrteexA z1 ;
qte11i = rteinA 0. z1 ;
qte12i = 0. rteinA z1 ;
qte13i = mrteinA 0. z1 ;
qte14i = 0. mrteinA z1 ;
 
**Points de l'interface 6 :
qte61e = rteexA 0. z6 ;
qte61i = rteinA 0. z6 ;
 
**Points de l'interface 7 :
qte71e = rteexG 0. z7 ;
qte72e = 0. rteexG z7 ;
qte73e = mrteexG 0. z7 ;
qte74e = 0. mrteexG z7 ;
qte71i = rteinG 0. z7 ;
qte72i = 0. rteinG z7 ;
qte73i = mrteinG 0. z7 ;
qte74i = 0. mrteinG z7 ;
 
* Points supplementaires pour la cavite de gaz :
qte21i = rteinA 0. z2 ;
qte31i = rteinA 0. z3 ;
qte41i = rteinA 0. z4 ;
qte51i = rteinA 0. z5 ;
qte81i = rteinG 0. z8 ;
qte91i = rteinG 0. z9 ;
qte101i = rteinG 0. z10 ;
qte111i = rteinG 0. z11 ;
qteM1i = rteinG 0. zM ;
 
*** Lignes des interfaces 1 -> 6 :
lte11e_p = 'CER3' np qte11e qte12e qte13e ;
lte12e_p = 'CER3' np qte13e qte14e qte11e ;
ltet1e_p = 'ET' lte11e_p lte12e_p ;
lte11i_p = 'CER3' np qte11i qte12i qte13i ;
lte12i_p = 'CER3' np qte13i qte14i qte11i ;
ltet1i_p = 'ET' lte11i_p lte12i_p ;
 
ltet2e_p ltet3e_p ltet4e_p ltet5e_p ltet6e_p no7 no8 no9 no10 no11 noM
 = PLUSJU ltet1e_p ;
ltet2i_p ltet3i_p ltet4i_p ltet5i_p ltet6i_p no7 no8 no9 no10 no11 noM
 = PLUSJU ltet1i_p ;
 
*** Lignes des interfaces 7 -> 11 :
lte71e_p = 'CER3' np qte71e qte72e qte73e ;
lte72e_p = 'CER3' np qte73e qte74e qte71e ;
ltet7e_p = 'ET' lte71e_p lte72e_p ;
lte71i_p = 'CER3' np qte71i qte72i qte73i ;
lte72i_p = 'CER3' np qte73i qte74i qte71i ;
ltet7i_p = 'ET' lte71i_p lte72i_p ;
 
ltet8e_p = PLUS ltet7e_p (0. 0. (z8 - z7)) ;
ltet8i_p = PLUS ltet7i_p (0. 0. (z8 - z7)) ;
ltet9e_p = PLUS ltet7e_p (0. 0. (z9 - z7)) ;
ltet9i_p = PLUS ltet7i_p (0. 0. (z9 - z7)) ;
ltet10e_ = PLUS ltet7e_p (0. 0. (z10 - z7)) ;
ltet10i_ = PLUS ltet7i_p (0. 0. (z10 - z7)) ;
ltet11e_ = PLUS ltet7e_p (0. 0. (z11 - z7)) ;
ltet11i_ = PLUS ltet7i_p (0. 0. (z11 - z7)) ;
 
*** Surfaces de toute les interfaces du tube de force externe
stfext1 = ( 'REGLER' ltet1i_p ne ltet1e_p ) ;
stfext2 = ( 'REGLER' ltet2i_p ne ltet2e_p ) ;
stfext3 = ( 'REGLER' ltet3i_p ne ltet3e_p ) ;
stfext4 = ( 'REGLER' ltet4i_p ne ltet4e_p ) ;
stfext5 = ( 'REGLER' ltet5i_p ne ltet5e_p ) ;
stfext6 = ( 'REGLER' ltet6i_p ne ltet6e_p ) ;
stfext7 = ( 'REGLER' ltet7i_p ne ltet7e_p ) ;
stfext8 = ( 'REGLER' ltet8i_p ne ltet8e_p ) ;
stfext9 = ( 'REGLER' ltet9i_p ne ltet9e_p ) ;
stfext10 = ( 'REGLER' ltet10i_ ne ltet10e_ ) ;
stfext11 = ( 'REGLER' ltet11i_ ne ltet11e_ ) ;
stfextM = PLUS stfext11 (0. 0. (zM - z11)) ;
 
elim 1.e-8 (qte21i ET qte31i ET qte41i ET qte51i ET qte61e ET qte61i ET
            qte81i ET qte91i ET qte101i ET qte111i ET qteM1i ET stfext2
            ET stfext3 ET stfext4 ET stfext5 ET stfext6 ET stfext8 ET
            stfext9 ET stfext10 ET stfext11 ET stfextM) ;
 
 
*      IV.9. Lame de gaz statique : ESSAI51 :
 
* Contrairement au autres maillages la lame de gaz se construit par
* rotation autour de l'axe plutot que par translation selon z.
 
ltiA = qti11e DROI qti21e nhautA ;
lteA = qte21i DROI qte11i nhautA ;
lsouA = qte11i DROI qti11e 1 ;
lsurA = qti21e DROI qte21i 1 ;
SgazA = DALL lsouA ltiA lsurA lteA 'PLAN' ;
 
ltiB = qti21e DROI qti31e nhautB ;
lteB = qte31i DROI qte21i nhautB ;
lsouB = qte21i DROI qti21e 1 ;
lsurB = qti31e DROI qte31i 1 ;
SgazB = DALL lsouB ltiB lsurB lteB 'PLAN' ;
 
ltiC = qti31e DROI qti41e nhautC ;
lteC = qte41i DROI qte31i nhautC ;
lsouC = qte31i DROI qti31e 1 ;
lsurC = qti41e DROI qte41i 1 ;
SgazC = DALL lsouC ltiC lsurC lteC 'PLAN' ;
 
ltiD = qti41e DROI qti51e nhautD ;
lteD = qte51i DROI qte41i nhautD ;
lsouD = qte41i DROI qti41e 1 ;
lsurD = qti51e DROI qte51i 1 ;
SgazD = DALL lsouD ltiD lsurD lteD 'PLAN' ;
 
ltiE = qti51e DROI qti61e nhautE ;
lteE = qte61i DROI qte51i nhautE ;
lsouE = qte51i DROI qti51e 1 ;
lsurE = qti61e DROI qte61i 1 ;
SgazE = DALL lsouE ltiE lsurE lteE 'PLAN' ;
 
ltiF = qti61e DROI qti71e nhautF ;
lteF = qte71i DROI qte61i nhautF ;
lsouF = qte61i DROI qti61e 1 ;
lsurF = qti71e DROI qte71i 1 ;
SgazF = DALL lsouF ltiF lsurF lteF 'PLAN' ;
 
ltiG = qti71e DROI qti81e nhautG ;
lteG = qte81i DROI qte71i nhautG ;
lsouG = qte71i DROI qti71e 1 ;
lsurG = qti81e DROI qte81i 1 ;
SgazG = DALL lsouG ltiG lsurG lteG 'PLAN' ;
 
ltiH = qti81e DROI qti91e nhautH ;
lteH = qte91i DROI qte81i nhautH ;
lsouH = qte81i DROI qti81e 1 ;
lsurH = qti91e DROI qte91i 1 ;
SgazH = DALL lsouH ltiH lsurH lteH 'PLAN' ;
 
ltiI = qti91e DROI qti101e nhautI ;
lteI = qte101i DROI qte91i nhautI ;
lsouI = qte91i DROI qti91e 1 ;
lsurI = qti101e DROI qte101i 1 ;
SgazI = DALL lsouI ltiI lsurI lteI 'PLAN' ;
 
ltiJ = qti101e DROI qti111e nhautJ ;
lteJ = qte111i DROI qte101i nhautJ ;
lsouJ = qte101i DROI qti101e 1 ;
lsurJ = qti111e DROI qte111i 1 ;
SgazJ = DALL lsouJ ltiJ lsurJ lteJ 'PLAN' ;
 
ltiK = qti111e DROI qtiM1e nhautK ;
lteK = qteM1i DROI qte111i nhautK ;
lsouK = qte111i DROI qti111e 1 ;
lsurK = qtiM1e DROI qteM1i 1 ;
SgazK = DALL lsouK ltiK lsurK lteK 'PLAN' ;
 
Sgaz = SgazA ET SgazB ET SgazC ET SgazD ET SgazE ET SgazF ET SgazG ET
       SgazH ET SgazI ET SgazJ ET SgazK ;
ELIM 1.e-8 Sgaz ;
 
 
** V. Construction des volumes :
 
*      V.1. Gaine - Fluide chaud - Panier - PE - Fluide froid - tdf interne :
 
*** definition des elements dans une hauteur elementaire
bahe = sthe1 ;
hahe = stheM ;
bahe_p = sthe_p1 ;
hahe_p = sthe_pM ;
tabfac = 'TABLE'   ; comm gaine        ;
tabfac_q = 'TABLE' ; comm panier       ;
tabele = 'TABLE'   ; comm fluide chaud ;
tabfac_p = 'TABLE' ; comm PE           ;
tabfac_t = 'TABLE' ; comm tdf interne  ;
tabele_t = 'TABLE' ; comm fluide froid ;
icpt = 1 ;
 
 
*** Ecriture des tables de correspondances :
 
* 1.PartA :
mtheA mthe_pA mtinA mtifA mtcbA mtin_pA mtif_pA mtcb_pA
mtin_peA mtif_peA mt_peA mtin_tA mtif_tA mt_tfiA
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe1 sthe_p1 lt1i lt1e stcb1 lt1i_p lt1e_p stcb_p1
             lpet1i lpet1e stpe1 ltit1i_p ltit1e_p stfint1
             z1 z2 nhautA
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 2.PartB :
mtheB mthe_pB mtinB mtifB mtcbB mtin_pB mtif_pB mtcb_pB
mtin_peB mtif_peB mt_peB mtin_tB mtif_tB mt_tfiB
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe2 sthe_p2 lt2i lt2e stcb2 lt2i_p lt2e_p stcb_p2
             lpet2i lpet2e stpe2 ltit2i_p ltit2e_p stfint2
             z2 z3 nhautB
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 3.PartC :
mtheC mthe_pC mtinC mtifC mtcbC mtin_pC mtif_pC mtcb_pC
mtin_peC mtif_peC mt_peC mtin_tC mtif_tC mt_tfiC
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe3 sthe_p3 lt3i lt3e stcb3 lt3i_p lt3e_p stcb_p3
             lpet3i lpet3e stpe3 ltit3i_p ltit3e_p stfint3
             z3 z4 nhautC
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 4.PartD :
mtheD mthe_pD mtinD mtifD mtcbD mtin_pD mtif_pD mtcb_pD
mtin_peD mtif_peD mt_peD mtin_tD mtif_tD mt_tfiD
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe4 sthe_p4 lt4i lt4e stcb4 lt4i_p lt4e_p stcb_p4
             lpet4i lpet4e stpe4 ltit4i_p ltit4e_p stfint4
             z4 z5 nhautD
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 5.PartE :
mtheE mthe_pE mtinE mtifE mtcbE mtin_pE mtif_pE mtcb_pE
mtin_peE mtif_peE mt_peE mtin_tE mtif_tE mt_tfiE
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe5 sthe_p5 lt5i lt5e stcb5 lt5i_p lt5e_p stcb_p5
             lpet5i lpet5e stpe5 ltit5i_p ltit5e_p stfint5
             z5 z6 nhautE
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 6.PartF :
mtheF mthe_pF mtinF mtifF mtcbF mtin_pF mtif_pF mtcb_pF
mtin_peF mtif_peF mt_peF mtin_tF mtif_tF mt_tfiF
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe6 sthe_p6 lt6i lt6e stcb6 lt6i_p lt6e_p stcb_p6
             lpet6i lpet6e stpe6 ltit6i_p ltit6e_p stfint6
             z6 z7 nhautF
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
 
* 7.PartG :
** Parties droites :
mtheG no1 mtinG mtifG mtcbG mtin_pG mtif_pG mtcb_pG
mtin_peG mtif_peG mt_peG no2 no3 no4
tabfac tabfac_q tabfac_p no5 tabele no6 icpt1
    = VOLUJU sthe7 sthe_p7 lt7i lt7e stcb7 lt7i_p lt7e_p stcb_p7
             lpet7i lpet7e stpe7 ltit7i_p ltit7e_p stfint7
             z7 z8 nhautG
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
** Changement de section et raccordement:
C1 = 0. 0. z7 ;
C2 = 0. 0. z8 ;
li = qti71i DROI 1 qti81i ;
le = qti71e DROI 1 qti81e ;
np2 = 2 * np ;
 
mthe_pG = sthe_p7 VOLU sthe_p8 1   ; comm remplace no1 ;
mtin_tG = li ROTA (np2) 360. C1 C2 ; comm remplace no2 ;
mtif_tG = le ROTA (np2) 360. C1 C2 ; comm remplace no3 ;
mt_tfiG = stfint7 VOLU stfint8 1   ; comm remplace no4 ;
 
nmtif1 = 'NBEL' mtin_tG ;
'REPETER' imtif1 nmtif1 ;
   lmtif1 = 'ELEM' mtin_tG &imtif1 ;
   tabfac_t . icpt = lmtif1 ;
   tabele_t . icpt = mthe_pG ;
   icpt = '+' icpt 1 ;
'FIN' imtif1 ;
 
 
* 8.PartH :
mtheH mthe_pH mtinH mtifH mtcbH mtin_pH mtif_pH mtcb_pH
mtin_peH mtif_peH mt_peH mtin_tH mtif_tH mt_tfiH
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe8 sthe_p8 lt8i lt8e stcb8 lt8i_p lt8e_p stcb_p8
             lpet8i lpet8e stpe8 ltit8i_p ltit8e_p stfint8
             z8 z9 nhautH
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 9.PartI : on ne modelise pas l'arret du panier et de la lame d'eau
*                              le changement de section du PE et des jambe
*                              le capteur ac=vec les bobines.
 
mtheI mthe_pI mtinI mtifI mtcbI mtin_pI mtif_pI mtcb_pI
mtin_peI mtif_peI mt_peI mtin_tI mtif_tI mt_tfiI
tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe9 sthe_p9 lt9i lt9e stcb9 lt9i_p lt9e_p stcb_p9
             lpet9i lpet9e stpe9 ltit9i_p ltit9e_p stfint9
             z9 z10 nhautI
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
* 10.PartJ :
** Parties droites :
no1 mthe_pJ no2 no3 no4 mtin_pJ mtif_pJ mtcb_pJ
mtin_peJ mtif_peJ mt_peJ mtin_tJ mtif_tJ mt_tfiJ
no5 tabfac_q tabfac_p tabfac_t no6 tabele_t icpt1
    = VOLUJU sthe10 sthe_p10 lt10i lt10e stcb10 lt10i_p lt10e_p stcb_p10
             lpet10i lpet10e stpe10 ltit10i_ ltit10e_ stfint10
             z10 z11 nhautJ
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
** Maillage et table du fluide.
mtheJ = sthe10 VOLU sthe11 1   ; comm remplace no1 ;
 
nmtif1 = 'NBEL' mtin_pJ ;
'REPETER' imtif1 nmtif1 ;
   tabele . icpt = mtheJ ; comm remplace no6 ;
   icpt = '+' icpt 1 ;
'FIN' imtif1 ;
 
* 11.PartK :
mtheK mthe_pK no2 no3 no4 mtin_pK mtif_pK mtcb_pK
mtin_peK mtif_peK mt_peK mtin_tK mtif_tK mt_tfiK
no5 tabfac_q tabfac_p tabfac_t tabele tabele_t icpt
    = VOLUJU sthe11 sthe_p11 lt11i lt11e stcb11 lt11i_p lt11e_p stcb_p11
             lpet11i lpet11e stpe11 ltit11i_ ltit11e_ stfint11
             z11 zM nhautK
             tabfac tabfac_q tabfac_p tabfac_t tabele tabele_t icpt ;
 
 
* 12. Assemblage :
 
mthe    = mtheA ET mtheB ET mtheC ET mtheD ET mtheE ET mtheF ET mtheG ET
          mtheH ET mtheI ET mtheJ ET mtheK ;
mthe_p  = mthe_pA ET mthe_pB ET mthe_pC ET mthe_pD ET mthe_pE ET mthe_pF
          ET mthe_pG ET mthe_pH ET mthe_pI ET mthe_pJ ET mthe_pK ;
mtin    = mtinA ET mtinB ET mtinC ET mtinD ET mtinE ET mtinF ET mtinG ET
          mtinH ET mtinI ;
mtif    = mtifA ET mtifB ET mtifC ET mtifD ET mtifE ET mtifF ET mtifG ET
          mtifH ET mtifI ;
mtcb    = mtcbA ET mtcbB ET mtcbC ET mtcbD ET mtcbE ET mtcbF ET mtcbG ET
          mtcbH ET mtcbI ;
mtin_p  = mtin_pA ET mtin_pB ET mtin_pC ET mtin_pD ET mtin_pE ET mtin_pF
          ET mtin_pG ET mtin_pH ET mtin_pI ET mtin_pJ ET mtin_pK ;
mtif_p  = mtif_pA ET mtif_pB ET mtif_pC ET mtif_pD ET mtif_pE ET mtif_pF
          ET mtif_pG ET mtif_pH ET mtif_pI ET mtif_pJ ET mtif_pK ;
mtcb_p  = mtcb_pA ET mtcb_pB ET mtcb_pC ET mtcb_pD ET mtcb_pE ET mtcb_pF
          ET mtcb_pG ET mtcb_pH ET mtcb_pI ET mtcb_pJ ET mtcb_pK ;
mtin_pe = mtin_peA ET mtin_peB ET mtin_peC ET mtin_peD ET mtin_peE ET
          mtin_peF ET mtin_peG ET mtin_peH ET mtin_peI ET mtin_peJ ET
          mtin_peK ;
mtif_pe = mtif_peA ET mtif_peB ET mtif_peC ET mtif_peD ET mtif_peE ET
          mtif_peF ET mtif_peG ET mtif_peH ET mtif_peI ET mtif_peJ ET
          mtif_peK ;
mt_pe   = mt_peA ET mt_peB ET mt_peC ET mt_peD ET mt_peE ET mt_peF ET
          mt_peG ET mt_peH ET mt_peI ET mt_peJ ET mt_peK ;
mtin_t  = mtin_tA ET mtin_tB ET mtin_tC ET mtin_tD ET mtin_tE ET mtin_tF
          ET mtin_tG ET mtin_tH ET mtin_tI ET mtin_tJ ET mtin_tK ;
mtif_t  = mtif_tA ET mtif_tB ET mtif_tC ET mtif_tD ET mtif_tE ET mtif_tF
          ET mtif_tG ET mtif_tH ET mtif_tI ET mtif_tJ ET mtif_tK ;
mt_tfi  = mt_tfiA ET mt_tfiB ET mt_tfiC ET mt_tfiD ET mt_tfiE ET mt_tfiF
          ET mt_tfiG ET mt_tfiH ET mt_tfiI ET mt_tfiJ ET mt_tfiK ;
 
ELIM 1.e-8 (mthe ET mthe_p ET mtin ET mtif ET mtcb ET mtin_p ET mtif_p
            ET mtcb_p ET mtin_pe ET mtif_pe ET mt_pe ET mtin_t ET
            mtif_t ET mt_tfi);
 
mtt = mthe ET mthe_p ET mtin ET mtif ET mtcb ET mtin_p ET mtif_p ET
      mtcb_p ET mtin_pe ET mtif_pe ET mt_pe ET mtin_t ET mtif_t ET
      mt_tfi ;
mtt2 = mtt AFFI (1./8) (0. 0. 0.) (0. 0. 1.) ;
 
 
*      V.2. Rallonges inferieure et superieure :
mtcbAr = VOLU stcb1r nhautA TRAN (0. 0. (z2 - z1));
 
mtcbErp = VOLU stcb5r nhautE TRAN (0. 0. (z6 - z5)) ;
mtcbFr  = VOLU stcb6r nhautF TRAN (0. 0. (z7 - z6)) ;
mtcbGr  = VOLU stcb7r nhautG TRAN (0. 0. (z8 - z7)) ;
mtcbHr  = VOLU stcb8r nhautH TRAN (0. 0. (z9 - z8)) ;
mtcbIr  = VOLU stcb9r nhautI TRAN (0. 0. (z10 - z9)) ;
mtcbEr = mtcbErp ET mtcbFr ET mtcbGr ET mtcbHr ET mtcbIr ;
 
ELIM 1.e-8 (mtcbEr ET stcb1r ET stcb2r ET stcb3r ET stcb4r ET stcb5r ET
            stcb6r ET stcb7r ET stcb8r ET stcb9r ET stcb10r) ;
 
 
*      V.3. Lame d'eau stagnante :
mt_lamA = VOLU st_lam1 nhautA TRAN (0. 0. (z2 - z1))    ;
mt_lamB = VOLU st_lam2 nhautB TRAN (0. 0. (z3 - z2))    ;
mt_lamC = VOLU st_lam3 nhautC TRAN (0. 0. (z4 - z3))    ;
mt_lamD = VOLU st_lam4 nhautD TRAN (0. 0. (z5 - z4))    ;
mt_lamE = VOLU st_lam5 nhautE TRAN (0. 0. (z6 - z5))    ;
mt_lamF = VOLU st_lam6 nhautF TRAN (0. 0. (z7 - z6))    ;
mt_lamG = VOLU st_lam7 nhautG TRAN (0. 0. (z8 - z7))    ;
mt_lamH = VOLU st_lam8 nhautH TRAN (0. 0. (z9 - z8))    ;
mt_lamI = VOLU st_lam9 nhautI TRAN (0. 0. (z10 - z9))   ;
mt_lamJ = VOLU st_lam10 nhautJ TRAN (0. 0. (z11 - z10)) ;
mt_lamK = VOLU st_lam11 nhautK TRAN (0. 0. (zM - z11))  ;
 
mt_lam = mt_lamA ET mt_lamB ET mt_lamC ET mt_lamD ET mt_lamE ET mt_lamF
         ET mt_lamG ET mt_lamH ET mt_lamI ET mt_lamJ ET mt_lamK ;
 
ELIM 1.e-8 (mt_lam ET st_lam1 ET st_lam2 ET st_lam3 ET st_lam4 ET
            st_lam5 ET st_lam6 ET st_lam7 ET st_lam8 ET st_lam9 ET
            st_lam10 ET st_lam11 ET st_lamM) ;
 
 
*      V.4. Tube de force externe :
 
*** Partie A->E :
mt_tfeA = VOLU stfext1 nhautA TRAN (0. 0. (z2 - z1)) ; comm volume ;
mt_tfeB = VOLU stfext2 nhautB TRAN (0. 0. (z3 - z2)) ;
mt_tfeC = VOLU stfext3 nhautC TRAN (0. 0. (z4 - z3)) ;
mt_tfeD = VOLU stfext4 nhautD TRAN (0. 0. (z5 - z4)) ;
mt_tfeE = VOLU stfext5 nhautE TRAN (0. 0. (z6 - z5)) ;
 
mten_tA = ltet1i_p TRAN nhautA (0. 0. (z2 - z1)) ; comm paroi interne ;
mten_tB = ltet2i_p TRAN nhautB (0. 0. (z3 - z2)) ;
mten_tC = ltet3i_p TRAN nhautC (0. 0. (z4 - z3)) ;
mten_tD = ltet4i_p TRAN nhautD (0. 0. (z5 - z4)) ;
mten_tE = ltet5i_p TRAN nhautE (0. 0. (z6 - z5)) ;
 
mtef_tA = ltet1e_p TRAN nhautA (0. 0. (z2 - z1)) ; comm paroi externe ;
mtef_tB = ltet2e_p TRAN nhautB (0. 0. (z3 - z2)) ;
mtef_tC = ltet3e_p TRAN nhautC (0. 0. (z4 - z3)) ;
mtef_tD = ltet4e_p TRAN nhautD (0. 0. (z5 - z4)) ;
mtef_tE = ltet5e_p TRAN nhautE (0. 0. (z6 - z5)) ;
 
*** Raccordement Partie F :
C21 = 0. 0. z6 ;
C22 = 0. 0. z7 ;
li2 = qte61i DROI nhautF qte71i ;
le2 = qte61e DROI nhautF qte71e ;
np2 = 2 * np ;
 
mt_tfeF = stfext6 VOLU stfext7 nhautF ;
mten_tF = li2 ROTA np2 360 C21 C22 ;
mtef_tF = le2 ROTA np2 360 C21 C22 ;
 
***Parties G->J :
mt_tfeG = VOLU stfext7 nhautG TRAN (0. 0. (z8 - z7)) ; comm volume ;
mt_tfeH = VOLU stfext8 nhautH TRAN (0. 0. (z9 - z8)) ;
mt_tfeI = VOLU stfext9 nhautI TRAN (0. 0. (z10 - z9)) ;
mt_tfeJ = VOLU stfext10 nhautJ TRAN (0. 0. (z11 - z10)) ;
mt_tfeK = VOLU stfext11 nhautK TRAN (0. 0. (zM - z11)) ;
 
mten_tG = ltet7i_p TRAN nhautG (0. 0. (z8 - z7)) ; comm paroi interne ;
mten_tH = ltet8i_p TRAN nhautH (0. 0. (z9 - z8)) ;
mten_tI = ltet9i_p TRAN nhautI (0. 0. (z10 - z9)) ;
mten_tJ = ltet10i_ TRAN nhautJ (0. 0. (z11 - z10)) ;
mten_tK = ltet11i_ TRAN nhautK (0. 0. (zM - z11)) ;
 
mtef_tG = ltet7e_p TRAN nhautG (0. 0. (z8 - z7)) ; comm paroi externe ;
mtef_tH = ltet8e_p TRAN nhautH (0. 0. (z9 - z8)) ;
mtef_tI = ltet9e_p TRAN nhautI (0. 0. (z10 - z9)) ;
mtef_tJ = ltet10e_ TRAN nhautJ (0. 0. (z11 - z10)) ;
mtef_tK = ltet11e_ TRAN nhautK (0. 0. (zM - z11)) ;
 
*** Assemblage :
mt_tfe = mt_tfeA ET mt_tfeB ET mt_tfeC ET mt_tfeD ET mt_tfeE ET mt_tfeF
          ET mt_tfeG ET mt_tfeH ET mt_tfeI ET mt_tfeJ ET mt_tfeK ;
mten_t = mten_tA ET mten_tB ET mten_tC ET mten_tD ET mten_tE ET mten_tF
          ET mten_tG ET mten_tH ET mten_tI ET mten_tJ ET mten_tK ;
mtef_t = mtef_tA ET mtef_tB ET mtef_tC ET mtef_tD ET mtef_tE ET mtef_tF
          ET mtef_tG ET mtef_tH ET mtef_tI ET mtef_tJ ET mtef_tK ;
 
ELIM 1.e-8 (mt_tfe et mtef_t et mten_t) ;
 
 
*      V.5. Lame de gaz :
 
*** V.5.1. Cavite de gaz :
* Surface interne de la lame de gaz = surface externe du tdf interne
* Orientation vers le gaz ie vers l'exterieur :
cavint = INVE (ORIE mtif_t 'POINT' (0. 0. zM)) ;
vers cavint ;
 
* Surface externe de la lame de gaz = surface interne du tdf externe
* Orientation vers le gaz ie vers l'interieur :
cavext = ORIE mten_t 'POINT' (0. 0. zM) ;
vers cavext ;
 
lcav31 = cavint et cavext ;
 
 
*** V.5.2. Volume interieur pour la conduction :
mt_gaz = Sgaz VOLU 'ROTA' (np2) 360. C1 C2 ;
 
 
 
 
*      V.6. ELIMINATION :
 
ELIM 1.e-8 (bahe ET hahe ET bahe_p ET hahe_p ET mthe ET mtin ET mtif ET
            mtcb ET mtin_p ET mtif_p ET mtcb_p ET mtin_pe ET mtif_pe ET
            mt_pe ET mtin_t ET mtif_t ET mt_tfi ET mtcbAr ET mtcbEr ET
            mt_lam ET mt_tfe ET mtef_t ET mten_t ET lcav31 ET mt_gaz) ;
 
**************************FIN DU MAILLAGE*******************************
************************************************************************
 
 
* CONSTRUCTION D'OBJETS UTILES AU CALCUL
************************************************************************
 
** I. Rayonnement
 
*      I.1. Description
* stefan  : constante de Stefan
* emistf  : emmissivite du tube de force
* emistfi : emmisivite du tube de force interne
* emistfe : emmisivite du tube de force externe
* alpha1  : emmisivite equivalente pour le face a face
 
*      I.2. Valeurs
stefan  = 5.67E-8 ;
emistf = 0.7 ;
emistfi = emistf ;
emistfe = emistf ;
r_1 = 1. - emistfi ;
r_2 = 1. - emistfe ;
alpha1 = (emistfi*emistfe) / (1.-(r_1*r_2)) ;
 
*      I.3. Rayonnement
mrcav1 = MODE lcav31 thermique RAYONNEMENT  'CAVITE';
ecav1  = MATE mrcav1 'EMIS' emistf ;
 
*      I.4 Parametres pour traiter le rayonnement.
*Nombre d'iterations internes et facteur de relaxation ;
SI rayo ;
 nitr1 = 10 ;
 omgr1 = 0.5 ;
SINON ;
 nitr1 = 1 ;
 omgr1 = 1. ;
FINSI ;
 
 
** II.Generation des QUAFs
 
_mtcbAr = 'CHANGER' mtcbAr 'QUAF' ;
_mtcbEr = 'CHANGER' mtcbEr 'QUAF' ;
_mtcb = 'CHANGER' mtcb 'QUAF' ;
_mtin = 'CHANGER' mtin 'QUAF' ;
_mtif = 'CHANGER' mtif 'QUAF' ;
 
_mthe = 'CHANGER' mthe 'QUAF' ;
_bahe = 'CHANGER' bahe 'QUAF' ;
_hahe = 'CHANGER' hahe 'QUAF' ;
 
_mtcb_p = 'CHANGER' mtcb_p 'QUAF' ;
_mtin_p = 'CHANGER' mtin_p 'QUAF' ;
_mtif_p = 'CHANGER' mtif_p 'QUAF' ;
_mt_lam = 'CHANGER' mt_lam 'QUAF' ;
_mt_pe = 'CHANGER' mt_pe 'QUAF' ;
_mtin_pe = 'CHANGER' mtin_pe 'QUAF' ;
_mtif_pe = 'CHANGER' mtif_pe 'QUAF' ;
 
_mthe_p = 'CHANGER' mthe_p 'QUAF' ;
_bahe_p = 'CHANGER' bahe_p 'QUAF' ;
_hahe_p = 'CHANGER' hahe_p 'QUAF' ;
 
_mt_tfi = 'CHANGER' mt_tfi 'QUAF' ;
_mtin_t = 'CHANGER' mtin_t 'QUAF' ;
_mtif_t = 'CHANGER' mtif_t 'QUAF' ;
_mt_tfe = 'CHANGER' mt_tfe 'QUAF' ;
_mten_t = 'CHANGER' mten_t 'QUAF' ;
_mtef_t = 'CHANGER' mtef_t 'QUAF' ;
_mt_gaz = 'CHANGER' mt_gaz 'QUAF' ;
 
_lcav31 = 'CHANGER' lcav31 'QUAF' ;
 
 
** III. CREATION DE LA MATRICE DE COUPLAGE
 
* Explications :
*    Generation de connectivites face-centre cfc et des matrices de
*    passage centre-face pour projeter T0 et face-centre pour ramener
*    une integrale de flux a une puissance volumique.
 
 
*      III.1. Interfaces Crayon - Fluide chaud
 
** Construction des matrices :
nitab = 'DIME' tabfac ;
lng = PROG ;
'REPETER' iitab  nitab ;
   ifac = tabfac . &iitab ;
   iqfac = 'CHANGER' ifac 'QUAF' ;
   miqfac = 'MODELISER' iqfac 'NAVIER_STOKES' DISCR ;
   ciqfac = ('DOMA' miqfac 'CENTRE') 'POIN' 1 ;
   iele = tabele . &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) ;
   lng = ( prog ( mesu clfc )) ;
        cfc = clfc ;
        mcf = mlcf ;
        mfc = mlfc ;
   'SINON' ;
   lng = lng et ( prog ( mesu clfc ) ) ;
      cfc = 'ET' cfc clfc ;
      mcf = 'ET' mcf mlcf ;
      mfc = 'ET' mfc mlfc ;
   'FINSI' ;
'FIN' iitab ;
 
** Elimination des points doubles :
mt = ('ET' mthe mtcb) ;
pelim = '*' ('**' ('/' ('MESURE' mt) ('NBEL' mt))
                  ('/' 1.D0 ('VALEUR' 'DIME')))
            1.D-3 ;
 
'ELIMINATION' (_bahe 'ET' _mthe 'ET' cfc) pelim ;
'ELIMINATION' (_mtcb 'ET' _mtin 'ET' _mtif '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') ;
 
 
*      III.2. Interfaces Fluide chaud - Panier
 
** Construction des matrices :
nitab = 'DIME' tabfac_q ;
lng = PROG ;
'REPETER' iitab  nitab ;
   ifac = tabfac_q . &iitab ;
   iqfac = 'CHANGER' ifac 'QUAF' ;
 
   miqfac = 'MODELISER' iqfac 'NAVIER_STOKES' DISCR ;
   ciqfac = ('DOMA' miqfac 'CENTRE') 'POIN' 1 ;
   iele = tabele . &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_q = 'MANUEL' 'SEG2' ciqfac ciqele ;
   mlcf_q = 'MANUEL' 'RIGIDITE' clfc_q ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 1.D0) ('PROG' 0.D0 0.D0) ;
   mlfc_q = 'MANUEL' 'RIGIDITE' clfc_q ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 0.D0)
            ('PROG' ('/' 1.D0 volele) 0.D0) ;
   'SI' ('EGA' &iitab 1) ;
   lng = ( prog ( mesu clfc_q )) ;
      cfc_q = clfc_q ;
      mcf_q = mlcf_q ;
      mfc_q = mlfc_q ;
   'SINON' ;
   lng = lng et ( prog ( mesu clfc_q ) ) ;
      cfc_q = 'ET' cfc_q clfc_q ;
      mcf_q = 'ET' mcf_q mlcf_q ;
      mfc_q = 'ET' mfc_q mlfc_q ;
   'FINSI' ;
 
'FIN' iitab ;
 
mess (mini lng) ( maxi lng ) ;
 
** Elimination de spoints doubles :
mt_q = ('ET' mthe mtcb_p) ;
pelim_q = '*' ('**' ('/' ('MESURE' mt_q) ('NBEL' mt_q))
                  ('/' 1.D0 ('VALEUR' 'DIME')))
            1.D-3 ;
 
'ELIMINATION' (_bahe 'ET' _mthe 'ET' cfc_q) pelim_q ;
'ELIMINATION' (_mtcb_p 'ET' _mtin_p 'ET' _mtif_p 'ET' cfc_q) pelim_q ;
mcf_q = 'CHANGER' mcf_q 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc_q = 'CHANGER' mfc_q 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
 
 
*     III.3. Interface PE - fluide froid
 
** Construction des matrices :
nitab = 'DIME' tabfac_p ;
lng = PROG ;
'REPETER' iitab  nitab ;
   ifac = tabfac_p . &iitab ;
   iqfac = 'CHANGER' ifac 'QUAF' ;
   miqfac = 'MODELISER' iqfac 'NAVIER_STOKES' DISCR ;
   ciqfac = ('DOMA' miqfac 'CENTRE') 'POIN' 1 ;
   iele = tabele_t . &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_p = 'MANUEL' 'SEG2' ciqfac ciqele ;
     mlcf_p = 'MANUEL' 'RIGIDITE' clfc_p ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 1.D0) ('PROG' 0.D0 0.D0) ;
     mlfc_p = 'MANUEL' 'RIGIDITE' clfc_p ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 0.D0)
            ('PROG' ('/' 1.D0 volele) 0.D0) ;
   'SI' ('EGA' &iitab 1) ;
   lng = ( prog ( mesu clfc )) ;
      cfc_p = clfc_p ;
      mcf_p = mlcf_p ;
      mfc_p = mlfc_p ;
   'SINON' ;
   lng = lng et ( prog ( mesu clfc_p ) ) ;
      cfc_p = 'ET' cfc_p clfc_p ;
      mcf_p = 'ET' mcf_p mlcf_p ;
      mfc_p = 'ET' mfc_p mlfc_p ;
   'FINSI' ;
'FIN' iitab ;
 
** Elimination des points doubles :
mt_p = ('ET' mthe_p mt_pe) ;
pelim_p = '*' ('**' ('/' ('MESURE' mt_p) ('NBEL' mt_p))
                  ('/' 1.D0 ('VALEUR' 'DIME')))
            1.D-3 ;
 
'ELIMINATION' (_bahe_p 'ET' _mthe_p 'ET' cfc_p) pelim_p ;
'ELIMINATION' (_mt_pe 'ET' _mtin_pe 'ET' _mtif_pe 'ET' cfc_p) pelim_p ;
mcf_p = 'CHANGER' mcf_p 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc_p = 'CHANGER' mfc_p 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
 
*     III.4. Interface tube force interne - fluide froid
 
** Construction des matrices :
nitab = 'DIME' tabfac_t ;
lng = PROG ;
'REPETER' iitab  nitab ;
   ifac = tabfac_t . &iitab ;
   iqfac = 'CHANGER' ifac 'QUAF' ;
   miqfac = 'MODELISER' iqfac 'NAVIER_STOKES' DISCR ;
   ciqfac = ('DOMA' miqfac 'CENTRE') 'POIN' 1 ;
   iele = tabele_t . &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_t = 'MANUEL' 'SEG2' ciqfac ciqele ;
   mlcf_t = 'MANUEL' 'RIGIDITE' clfc_t ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 1.D0) ('PROG' 0.D0 0.D0) ;
   mlfc_t = 'MANUEL' 'RIGIDITE' clfc_t ('MOTS' 'T') 'QUEL'
            ('PROG' 0.D0 0.D0)
            ('PROG' ('/' 1.D0 volele) 0.D0) ;
   'SI' ('EGA' &iitab 1) ;
   lng = ( prog ( mesu clfc_t )) ;
      cfc_t = clfc_t ;
      mcf_t = mlcf_t ;
      mfc_t = mlfc_t ;
   'SINON' ;
   lng = lng et ( prog ( mesu clfc_t ) ) ;
      cfc_t = 'ET' cfc_t clfc_t ;
      mcf_t = 'ET' mcf_t mlcf_t ;
      mfc_t = 'ET' mfc_t mlfc_t ;
   'FINSI' ;
'FIN' iitab ;
 
** Elimination des points doubles :
mt_t = ('ET' mthe_p mt_tfi) ;
pelim_t = '*' ('**' ('/' ('MESURE' mt_t) ('NBEL' mt_t))
                  ('/' 1.D0 ('VALEUR' 'DIME')))
            1.D-3 ;
 
'ELIMINATION' (_bahe_p 'ET' _mthe_p 'ET' cfc_t) pelim_t ;
'ELIMINATION' (_mt_tfi 'ET' _mtin_t 'ET' _mtif_t 'ET' cfc_t) pelim_t ;
mcf_t = 'CHANGER' mcf_t 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc_t = 'CHANGER' mfc_t 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
 
 
tcpu9 ;
 
********FIN DE LA CONSTRUCTION D'OBJETS UTILES AU CALCUL ***************
************************************************************************
 
 si interact ; opti donn 5 ; finsi ;
 
si rayo ;
*
* Elements de liaison de la cavite pour le rayonnement face a face
* ----------------------------------------------------------------
nepbth = 1 ; comm nombre elements radiaux de mt_gaz ;
rbz0 = rayon0 mt_gaz ;
epcav = (maxi rbz0 ) - ( mini rbz0 ) ;
 
opti dime 3 elem qua4 ;
repe x ( nbel cavint ) ;
 elep = cavint elem &x ;
 elep1 = chan poi1 elep ;
 pp1 = point 1 elep1;
 pgg1 = point cavext proc pp1 ;
 pgg2 = point cavext proc (point 2 elep1 ) ;
 pgg3 = point cavext proc (point 3 elep1 ) ;
 pgg4 = point cavext proc (point 4 elep1 ) ;
 ll1 = pgg1 d 1 pgg2 ;
 ll2 = pgg2 d 1 pgg3 ;
 ll3 = pgg3 d 1 pgg4 ;
 ll4 = pgg4 d 1 pgg1 ;
 elig1 = daller ll1 ll2 ll3 ll4 ;
 nouv0 = liai epcav elep elig1  ;
  new   = pp1 d 1 pgg1 ;
 si ( ega &x 1 ) ;
  rel12q = new   ;
  rac12 = nouv0 ;
 sinon ;
  rel12q = rel12q et new ;
  rac12 = rac12 et nouv0 ;
 finsi ;
fin x ;
mess (nbel cavext ) ( nbel cavint ) (nbel rac12 ) (nbel rel12q) ;
 
finsi ; comm si rayo ;
graph     = graph0 ;
graffic   = graffic0 ;
interact  = interac0 ;
rayo      = rayo0 ;
egasol = egasol0 ;
egaliq = egaliq0 ;
 
verbose   = verbose0 ;
precis    = precis0 ;
complet   = comple0 ;
 
prhe   = prhe0 ;
qmhe   = qmhe0 ;
theec  = theec0 ;
tpisci = tpisci0 ;
plin1  = plin10 ;
psatb  = psatb0 ;
ppa    = ppa0 ;
posix  = posix0 ;
Egmax = Egmax0 ;
 
varcna   = varcna0 ;
hcncst   = hcncst0 ;
 
si (NON complet ) ;
mess 'Complet = FAUX' ;
 nitr1 = 1 ;
 omgr1 = 1. ;
 nb_bouc = 1 ;
sinon ;
mess 'Complet = VRAI' ;
FINSI ;
 
theek   = theec + CK ;
tpiscik = tpisci + CK ;
plin1   = '*' plin1 100.D0 ;
Psat = PsatB * 1.e5 ;
Tsat = CK + ( h2o_tsat Psat );
 
************************************************************************
* Calcul des proprietes thermodynamiques de l'eau (procedure PRPEAU)   *
************************************************************************
opti echo 0 ;
prt =  ('PROG' CK 'PAS' 10. (tsat - 1. )) 'ET' ('PROG' 1000. ) ;
ro_o = prog ; mu_o = prog ; lamb_o = prog ; cp_o = prog ;
repe x ((dime prt ) - 1 ) ;
 ro_o   = ro_o   et ( prog  (h2o_rhol ((extr prt &x)- CK ) ppa )) ;
 mu_o   = mu_o   et ( prog  (h2o_mu9  ((extr prt &x)- CK ) ppa )) ;
 cp_o   = cp_o   et ( prog  (h2o_cplq ((extr prt &x)- CK ) ppa )) ;
 lamb_o = lamb_o et ( prog  (h2o_lmda ((extr prt &x)- CK ) ppa )) ;
fin x ;
* On  duplique le dernier element
* -------------------------------
ro_o = ro_o et (prog ( extr ro_o (dime ro_o ))) ;
mu_o = mu_o et (prog ( extr mu_o (dime mu_o ))) ;
cp_o = cp_o et (prog ( extr cp_o (dime cp_o ))) ;
lamb_o = lamb_o et (prog ( extr lamb_o (dime lamb_o ))) ;
list ro_o; list mu_o; list cp_o;list lamb_o ;
opti echo 1 ;
 
 
************************************************************************
************************************************************************
 
********** RAJOUT ESSAI14 :Analyse monophasique ************************
 
* PREPARATION DU CALCUL MONOPHASIQUE
************************************************************************
 
** I. Modeles Navier-Stokes
 
$mtcbAr = 'MODELISER' _mtcbAr 'NAVIER_STOKES' discr ;
$mtcbEr = 'MODELISER' _mtcbEr 'NAVIER_STOKES' discr ;
$mtcb = 'MODELISER' _mtcb 'NAVIER_STOKES' discr ;
$mtin = 'MODELISER' _mtin 'NAVIER_STOKES' discr ;
$mtif = 'MODELISER' _mtif 'NAVIER_STOKES' discr ;
 
$bahe = 'MODELISER' _bahe 'NAVIER_STOKES' discr ;
$hahe = 'MODELISER' _hahe 'NAVIER_STOKES' discr ;
$mthe = 'MODELISER' _mthe 'NAVIER_STOKES' discr ;
 
$mtcb_p = 'MODELISER' _mtcb_p 'NAVIER_STOKES' discr ;
$mtin_p = 'MODELISER' _mtin_p 'NAVIER_STOKES' discr ;
$mtif_p = 'MODELISER' _mtif_p 'NAVIER_STOKES' discr ;
$mt_lam = 'MODELISER' _mt_lam 'NAVIER_STOKES' discr ;
$mt_pe  = 'MODELISER' _mt_pe 'NAVIER_STOKES' discr ;
$mtin_pe = 'MODELISER' _mtin_pe 'NAVIER_STOKES' discr ;
$mtif_pe = 'MODELISER' _mtif_pe 'NAVIER_STOKES' discr ;
 
$bahe_p = 'MODELISER' _bahe_p 'NAVIER_STOKES' discr ;
$hahe_p = 'MODELISER' _hahe_p 'NAVIER_STOKES' discr ;
$mthe_p = 'MODELISER' _mthe_p 'NAVIER_STOKES' discr ;
 
$mt_tfi  = 'MODELISER'  _mt_tfi 'NAVIER_STOKES' discr ;
$mt_tfe  = 'MODELISER'  _mt_tfe 'NAVIER_STOKES' discr ;
$mt_gaz  = 'MODELISER'  _mt_gaz 'NAVIER_STOKES' discr ;
$mtin_t  = 'MODELISER'  _mtin_t 'NAVIER_STOKES' discr ;
$mten_t  = 'MODELISER'  _mten_t 'NAVIER_STOKES' discr ;
$mtif_t  = 'MODELISER'  _mtif_t 'NAVIER_STOKES' discr ;
$mtef_t  = 'MODELISER'  _mtef_t 'NAVIER_STOKES' discr ;
 
 si rayo ;
** II. Special rayonnement face a face
 
dcav_1   = 'DOMA' $mtif_t  'MAILLAGE' ;
dcav_2   = 'DOMA' $mten_t  'MAILLAGE' ;
mess (nbel rac12) (nbel dcav_1) ( nbel dcav_2 ) ;
 
*mess (nbno dcav_1) ( nbno dcav_2 ) ( nbno rac12 );
sr_12  = rac12  'CHANGER' 'LIGNE' ;
sr_1   = dcav_1 'CHANGER' 'LIGNE' ;
sr_2   = dcav_2 'CHANGER' 'LIGNE' ;
rel_12 = diff sr_12 ( sr_2 ET sr_1 ) ;
mess (nbno sr_12) (nbno sr_1)(nbno sr_2) ;
*trac ( sr_1 et (sr_2 coul roug) et (rel_12 coul vert));
 
finsi ; comm si rayo ;
 
** III. Calcul des parametres hydrauliques
 
**    III.1. Coefficients permettant d'affecter les valeurs des
*           grandeurs aux 3 parties differentes :
 
** III.1.1. Pour le fluide chaud :
dmthe = DOMA $mthe 'CENTRE' ;
zhe = COOR 3 dmthe ;
cA_c = zhe MASQ 'INFERIEUR' z10 ;
cK_c = zhe MASQ 'SUPERIEUR' z11 ;
cJ1_c = zhe MASQ 'SUPERIEUR' z10 ;
cJ2_c = zhe MASQ 'INFERIEUR' z11 ;
cJ_c = cJ1_c + cJ2_c - 1. ;
 
*Pour la paroi interne du panier interne :
dmtin_p = DOMA $mtin_p 'CENTRE' ;
zin_p = COOR 3 dmtin_p ;
 
cA_p = zin_p MASQ 'INFERIEUR' z10 ;
cK_p = zin_p MASQ 'SUPERIEUR' z11 ;
cJ1_p = zin_p MASQ 'SUPERIEUR' z10 ;
cJ2_p = zin_p MASQ 'INFERIEUR' z11 ;
cJ_p = cJ1_p + cJ2_p - 1. ;
 
 
** III.1.2. Pour le fluide froid :
dmthe_p = DOMA $mthe_p 'CENTRE' ;
zhe_p = COOR 3 dmthe_p ;
 
cA_f = zhe_p MASQ 'INFERIEUR' z7 ;
cH_f = zhe_p MASQ 'SUPERIEUR' z8 ;
cG1_f = zhe_p MASQ 'SUPERIEUR' z7 ;
cG2_f = zhe_p MASQ 'INFERIEUR' z8 ;
cG_f = cG1_f + cG2_f - 1. ;
 
*Pour la paroi externe du PE :
dmtif_pe = DOMA $mtif_pe 'CENTRE' ;
zif_pe = COOR 3 dmtif_pe ;
 
cA_pe = zif_pe MASQ 'INFERIEUR' z7 ;
cH_pe = zif_pe MASQ 'SUPERIEUR' z8 ;
cG1_pe = zif_pe MASQ 'SUPERIEUR' z7 ;
cG2_pe = zif_pe MASQ 'INFERIEUR' z8 ;
cG_pe = cG1_pe + cG2_pe - 1. ;
 
*Pour la paroi interne du tdf interne :
dmtin_t = DOMA $mtin_t 'CENTRE' ;
zin_t = COOR 3 dmtin_t ;
 
cA_t = zin_t MASQ 'INFERIEUR' z7 ;
cH_t = zin_t MASQ 'SUPERIEUR' z8 ;
cG1_t = zin_t MASQ 'SUPERIEUR' z7 ;
cG2_t = zin_t MASQ 'INFERIEUR' z8 ;
cG_t = cG1_t + cG2_t - 1. ;
 
*      III.2. Diametres hydrauliques exact-discretise
*Le diametre hydraulique se calcule, a partir de la formule :
*       " Dh = 4 * Spassage / Pmouille "
*SpassCf est une moyenne quadratique
 
** Fluide chaud : Dh est variable
*Parties A->I :
SpassAc = MAXI (RESU (DOMA $bahe 'VOLUME')) ;
Pif = (MAXI (RESU (DOMA $mtif 'VOLUME'))) / (z10 - z1) ;
Pin_p = KOPS (DOMA $mtin_p 'VOLUME') * cA_p ;
Pin_p = (MAXI (RESU Pin_p)) / (z10 - z1) ;
diamhAc = 4 * SpassAc / (Pif + Pin_p) ;
 
*Partie K :
SpassKc = MAXI (RESU (DOMA $hahe 'VOLUME')) ;
Pin_p = KOPS (DOMA $mtin_p 'VOLUME') * cK_p ;
Pin_p = (MAXI (RESU Pin_p)) / (zM - z11) ;
diamhKc = 4 * SpassKc / Pin_p ;
 
*Partie J :
SpassJc = (((SpassAc ** 0.5) + (SpassKc ** 0.5)) / 2.) ** 2. ;
Pin_p = KOPS (DOMA $mtin_p 'VOLUME') * cJ_p ;
Pin_p = (MAXI (RESU Pin_p)) / (z11 - z10) ;
diamhJc = 4 * SpassJc / Pin_p ;
 
diamh = diamhAc ;
diamh_p = (KOPS cA_p * diamhAc) + (KOPS cK_p * diamhKc) +
           (KOPS cJ_p * diamhJc) ;
 
** Fluide froid : Dh est variable
*Parties A->F :
SpassA = MAXI (RESU (DOMA $bahe_p 'VOLUME')) ;
Pif_pe = KOPS (DOMA $mtif_pe 'VOLUME') * cA_pe ;
Pif_pe = (MAXI (RESU Pif_pe)) / (z7 - z1) ;
Pin_t = KOPS (DOMA $mtin_t 'VOLUME') * cA_t ;
Pin_t = (MAXI (RESU Pin_t)) / (z7 - z1) ;
diamhA = 4 * SpassA / (Pif_pe + Pin_t) ;
 
*Partie H :
SpassH = MAXI (RESU (DOMA $hahe_p 'VOLUME')) ;
Pif_pe = KOPS (DOMA $mtif_pe 'VOLUME') * cH_pe ;
Pif_pe = (MAXI (RESU Pif_pe)) / (zM - z8) ;
Pin_t = KOPS (DOMA $mtin_t 'VOLUME') * cH_t ;
Pin_t = (MAXI (RESU Pin_t)) / (zM - z8) ;
diamhH = 4 * SpassH / (Pif_pe + Pin_t) ;
 
*Partie G :
SpassG = (((SpassA ** 0.5) + (SpassH ** 0.5)) / 2.) ** 2. ;
Pif_pe = KOPS (DOMA $mtif_pe 'VOLUME') * cG_pe ;
Pif_pe = (MAXI (RESU Pif_pe)) / (z8 - z7) ;
Pin_t = KOPS (DOMA $mtin_t 'VOLUME') * cG_t ;
Pin_t = (MAXI (RESU Pin_t)) / (z8 - z7) ;
diamhG = 4 * SpassG / (Pif_pe + Pin_t) ;
 
diamh_pe = (KOPS cA_pe * diamhA) + (KOPS cH_pe * diamhH) +
           (KOPS cG_pe * diamhG) ;
diamh_t = (KOPS cA_t * diamhA) + (KOPS cH_t * diamhH) +
           (KOPS cG_t * diamhG) ;
 
 
*      III.3. Puissance surfacique dans le crayon chauffant
 
ptot = plin1 * (z4 - z3) ;
surfin2 = PI *  rint * 2.D0 * (z4 - z3)  ;
surfin =(MAXI (RESU (DOMA $mtin 'VOLUME'))) * ((z4 - z3) / (z10 - z1));
*'SI' ('NEG' surfin surfin2 (surfin * 1.D-2 )) ;
*   'MESSAGE' ('CHAINE' 'surfin=' surfin) ;
*   'MESSAGE' ('CHAINE' 'surfin2=' surfin2) ;
*   'ERREUR' 'Discretisation pas assez fine ?' ;
*'FINSI' ;
pmsurf = '/' ptot surfin ;
'MESSAGE' ('CHAINE' 'surfin  =' surfin ' m2') ;
'MESSAGE' ('CHAINE' 'surfin2 =' surfin2 ' m2') ;
'MESSAGE' ('CHAINE' 'pmsurf=' pmsurf ' W.m-2') ;
'MESSAGE' ('CHAINE' 'ptot=' ptot ' W') ;
 
* Coefficient permettant d'annuler la puissance hors du crayon
mprof = 'DOMA' $mtin 'CENTRE' ;
zcb = COOR 3 mprof ;
cW3 = zcb MASQUE 'SUPERIEUR' z3 ;
cW4 = zcb MASQUE 'INFERIEUR' z4 ;
cW  = KOPS (cW3 + cW4) - 1. ;
pmsurf = cW * pmsurf ;
 
vprof = COEFLU mprof ;
pmsurf = '*' vprof pmsurf ;
 
 
*      III.4. Debit des fluides chaud et froid.
* La vitesse massique est consideree constante par portion dans les fluides :
* AVANT : qhe = qmhe / Spass ;
qheAc = qmhe / SpassAc ;
qheJc = qmhe / SpassJc ;
qheKc = qmhe / SpassKc ;
qheA = qmhe / SpassA ;
qheG = qmhe / SpassG ;
qheH = qmhe / SpassH ;
 
qhe = (KOPS cA_c * qheAc) + (KOPS cJ_c * qheJc) + (KOPS cK_c * qheKc) ;
qhe_ga = qheAc ;
qhe_pa =(KOPS cA_p * qheAc) + (KOPS cJ_p * qheJc) + (KOPS cK_c * qheKc);
 
qhe_p = (KOPS cA_f * qheA) + (KOPS cG_f * qheG) + (KOPS cH_f * qheH) ;
qhe_pe =(KOPS cA_pe * qheA) + (KOPS cG_pe * qheG) + (KOPS cH_pe * qheH);
qhe_t = (KOPS cA_t * qheA) + (KOPS cG_t * qheG) + (KOPS cH_t * qheH) ;
 
 
*      III.5. Sources internes reparties  dans structure et/ou fluides
*    III.5.1. Description
*
* rhoeau0  : masse volumique de l'eau a la temperature d'entree
* rhozir0  : masse volumique du Zircaloy4
* rhohel0  : masse volumique de l'helium
* Egai      :| sources internes reparties  donne par la radioprotection (W/kg).
*           | = 0 si pas de puissance de sources internes reparties  dans l'eau.
* Wgaeau_e : puissance de puissance de sources internes de l'eau en entree (W/m3).
* Wgazir_e  : puissance de puissance de sources internes initiale du zircaloy (W/m3).
* Wgahel_e  : puissance de puissance de sources internes initiale de l'helium (W/m3)
*
*RQ : il faudra prendre des vrais tables pour les masses volumiques et
*     les echauffements
 
*    III.5.2. Valeurs :
 
rhoeau0 ju1 ju2 ju3 = PRPEAU prt ro_o mu_o cp_o lamb_o theek ;
rhozir0 ju1 = PRPZIRC theek ;
rhohel0 ju1 = PRPN2 theek ;
 
Egacb   = GAMMA zir posix $mtcb ;
Egacb_p = GAMMA zir posix $mtcb_p ;
Ega_tfi = GAMMA zir posix $mt_tfi ;
Ega_gaz = GAMMA hel posix $mt_gaz ;
Ega_tfe = GAMMA zir posix $mt_tfe ;
Egahe   = GAMMA eau posix $mthe ;
Egahe_p = GAMMA eau posix $mthe_p ;
 
**Rajout essai36 :
EgacbAr = GAMMA zir posix $mtcbAr ;
EgacbEr = GAMMA zir posix $mtcbEr ;
Ega_lam = GAMMA eau posix $mt_lam ;
Ega_pe = GAMMA zir posix $mt_pe ; comm a changer en inox ;
**
 
Wgacb0   = rhozir0 * Egacb ;
Wgacb_p0 = rhozir0 * Egacb_p ;
Wga_tfi0 = rhozir0 * Ega_tfi ;
Wga_gaz0 = rhohel0 * Ega_gaz ;
Wga_tfe0 = rhozir0 * Ega_tfe ;
Wgahe0   = rhoeau0 * Egahe ;
Wgahe_p0 = rhoeau0 * Egahe_p ;
 
**Rajout essai36 :
WgacbAr0 = rhozir0 * EgacbAr ;
WgacbEr0 = rhozir0 * EgacbEr ;
Wga_lam0 = rhoeau0 * Ega_lam ;
Wga_pe0 = rhozir0 * Ega_pe ; comm a changer en inox ;
**
 
** IV. Initialisation des tables pour les calculs thermiques
*      dans le crayon chauffant et dans l'eau
 
*      IV.1. Dans le crayon chauffant
 
rvcb = 'EQEX' 'ITMA' 1 'NITER' nitr1 'FIDT' 1000 'OMEGA' omgr1 ;
 
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
      'ZONE' $mtcbAr 'OPER' 'LAPN' 'LCBAr' 'INCO' 'TCB'
      'ZONE' $mtcbEr 'OPER' 'LAPN' 'LCBEr' 'INCO' 'TCB'
      'ZONE' $mtcb 'OPER' 'LAPN' 'LCB' 'INCO' 'TCB'
      'ZONE' $mtcb_p 'OPER' 'LAPN' 'LCB_p' 'INCO' 'TCB'
      'ZONE' $mt_lam 'OPER' 'LAPN' 'LCB_lam' 'INCO' 'TCB'
      'ZONE' $mt_pe 'OPER' 'LAPN' 'LCB_pe' 'INCO' 'TCB'
         ;
 
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
      'ZONE' $mt_tfi 'OPER' 'LAPN' 'LCB_tfi' 'INCO' 'TCB'
      'ZONE' $mt_tfe 'OPER' 'LAPN' 'LCB_tfe' 'INCO' 'TCB'
      'ZONE' $mt_gaz 'OPER' 'LAPN' 'LCB_gaz' 'INCO' 'TCB'
         ;
 
* Condition de convection naturelle a la piscine :
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtef_t 'OPER' 'ECHI' 'Hn4' 'Tp4' 'INCO' 'TCB' 'TCB'
          ;
 
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtin 'OPER' 'FIMP' 'PCB' 'INCO' 'TCB'
         ;
 
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtif 'OPER' 'ECHI' 'H' 'T0' 'INCO' 'TCB' 'TCB'
         'ZONE' $mtin_p 'OPER' 'ECHI' 'H1' 'T01' 'INCO' 'TCB' 'TCB'
         'ZONE' $mtif_pe 'OPER' 'ECHI' 'H2' 'T02' 'INCO' 'TCB' 'TCB'
         'ZONE' $mtin_t 'OPER' 'ECHI' 'H3' 'T03' 'INCO' 'TCB' 'TCB'
           ;
 
si egasol ;
***Puissance de sources internes dans les solides + lame de gaz
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcbAr 'OPER' 'FIMP' 'WgacbAr' 'INCO' 'TCB'
         'ZONE' $mtcbEr 'OPER' 'FIMP' 'WgacbEr' 'INCO' 'TCB'
         'ZONE' $mtcb 'OPER' 'FIMP' 'Wgacb' 'INCO' 'TCB'
         'ZONE' $mtcb_p 'OPER' 'FIMP' 'Wgacb_p' 'INCO' 'TCB'
         'ZONE' $mt_lam 'OPER' 'FIMP' 'Wga_lam' 'INCO' 'TCB'
         'ZONE' $mt_pe 'OPER' 'FIMP' 'Wga_pe' 'INCO' 'TCB'
         ;
 
rvcb = 'EQEX' rvcb
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mt_tfi 'OPER' 'FIMP' 'Wga_tfi' 'INCO' 'TCB'
         'ZONE' $mt_gaz 'OPER' 'FIMP' 'Wga_gaz' 'INCO' 'TCB'
         'ZONE' $mt_tfe 'OPER' 'FIMP' 'Wga_tfe' 'INCO' 'TCB'
         ;
finsi ;
 
***Rayonnement face a face
SI  rayo ;
rvcb = 'EQEX' rvcb
      'OPTI' 'EF' 'IMPL'
      'ZONE' $mtif_t 'OPER' 'ECHI' 'HR1' 'T2I' 'INCO' 'TCB'
      'ZONE' $mten_t 'OPER' 'ECHI' 'HR2' 'T1S' 'INCO' 'TCB'
        ;
FINSI ; comm rayonnement face a face ;
 
***Iterative solver for the fuel+structure part
*
* 1 : direct solver, more robust, slower
* 3 : iterative solver, less robust, faster
rvcb . 'METHINV' . 'TYPINV' = 1 ;
*rvcb . 'METHINV' . 'TYPINV' = 3 ;
 
*      IV.2. Dans l'eau
 
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'
'ZONE' $mthe_p 'OPER' 'KONV' 'RCPHE_p' 'UHE_p' 'LHE_p' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'LAPN' 'LHE' 'INCO' 'THE'
         'ZONE' $mthe_p 'OPER' 'LAPN' 'LHE_p' 'INCO' 'THE'
          ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'FIMP' 'PHE' 'INCO' 'THE'
         'ZONE' $mthe 'OPER' 'FIMP' 'PHE1' 'INCO' 'THE'
         'ZONE' $mthe_p 'OPER' 'FIMP' 'PHE2' 'INCO' 'THE'
         'ZONE' $mthe_p 'OPER' 'FIMP' 'PHE3' 'INCO' 'THE'
         ;
 
si egaliq ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'FIMP' 'Wgahe' 'INCO' 'THE'
         'ZONE' $mthe_p 'OPER' 'FIMP' 'Wgahe_p' 'INCO' 'THE'
         ;
finsi ;
 
***Conditions aux limites en temperature des jambes chaudes et froides
*   Ici on impose les temperatures en bas du crayon chauffant:
*       - jambe chaude: $bahe, Tch = theek + 1.
*       - jambe froide: $bahe_p, Tfr = theek
rvhe = 'EQEX' rvhe
         'CLIM' 'THE' 'TIMP' ('DOMA' $bahe 'MAILLAGE') theek
       'CLIM' 'THE' 'TIMP' ('DOMA' $bahe_p 'MAILLAGE') (theek - 1.)
          ;
*      IV.3. Initialisation de toutes les variables avant le calcul
 
rvi = 'TABLE' 'INCO' ;
rvcb . 'INCO' = rvi ;
rvhe . 'INCO' = rvi ;
rvi . 'TCB' = 'KCHT'
    ( $mtcbAr et $mtcbEr et $mtcb et $mtcb_p et $mt_lam et $mt_pe et
      $mt_tfi et $mt_tfe et $mt_gaz ) 'SCAL' 'SOMMET' theek ;
 
ju641 lcb = PRPZIRC theek ;
ju1 lcb_g = PRPN2 theek ; comm diazote ;
**Inox a rajouter
 
rvi . 'LCBAr' = 'KCHT' $mtcbAr 'SCAL' 'CENTRE' lcb ;
rvi . 'LCBEr' = 'KCHT' $mtcbEr 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB' = 'KCHT' $mtcb 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_p' = 'KCHT' $mtcb_p 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_lam' = 'KCHT' $mt_lam 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_pe' = 'KCHT' $mt_pe 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_tfi' = 'KCHT' $mt_tfi 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_tfe' = 'KCHT' $mt_tfe 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_gaz' = 'KCHT' $mt_gaz 'SCAL' 'CENTRE' lcb_g ;
 
rvi . 'PCB' = 'KCHT' $mtin 'SCAL' 'CENTRE' pmsurf ;
rvi . 'T0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T01' = 'KCHT'  $mtin_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T02' = 'KCHT'  $mtif_pe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T03' = 'KCHT'  $mtin_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F01' = 'KCHT'  $mtin_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F02' = 'KCHT'  $mtif_pe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F03' = 'KCHT'  $mtin_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H'  = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H1'  = 'KCHT'  $mtin_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H2'  = 'KCHT'  $mtif_pe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H3'  = 'KCHT'  $mtin_t 'SCAL' 'CENTRE' 0.D0 ;
SI rayo ;
rvi . 'T1S' = 'KCHT'  $mten_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T2I' = 'KCHT'  $mtif_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'HR2' = 'KCHT'  $mten_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'HR1' = 'KCHT'  $mtif_t 'SCAL' 'CENTRE' 0.D0 ;
FINSI ; comm rayonnement face a face ;
rvi . 'THE' = 'KCHT' ( $mthe ET $mthe_p ) 'SCAL' 'SOMMET' theek ;
rvi . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PHE1' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PHE2' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PHE3' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' 0.D0 ;
 
*Essai27 : changement de la CL a la piscine :
*------------------------------------------------
rvi . 'Hn4'  = 'KCHT' $mtef_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'Tp4' = 'KCHT' $mtef_t 'SCAL' 'CENTRE' tpiscik ;
rvi . 'Fn4' = 'KCHT' $mtef_t 'SCAL' 'CENTRE' 0.D0 ;
*------------------------------------------------
 
***Initialisation de la puissance des sources internes :
si egasol ;
 rvi . 'WgacbAr' = 'KCHT' $mtcbAr 'SCAL' 'CENTRE' WgacbAr0 ;
 rvi . 'WgacbEr' = 'KCHT' $mtcbEr 'SCAL' 'CENTRE' WgacbEr0 ;
 rvi . 'Wgacb' = 'KCHT' $mtcb 'SCAL' 'CENTRE' Wgacb0 ;
 rvi . 'Wgacb_p' = 'KCHT' $mtcb_p 'SCAL' 'CENTRE' Wgacb_p0 ;
 rvi . 'Wga_lam' = 'KCHT' $mt_lam 'SCAL' 'CENTRE' Wga_lam0 ;
 rvi . 'Wga_pe' = 'KCHT' $mt_pe 'SCAL' 'CENTRE' Wga_pe0 ;
 rvi . 'Wga_tfi' = 'KCHT' $mt_tfi 'SCAL' 'CENTRE' Wga_tfi0 ;
 rvi . 'Wga_gaz' = 'KCHT' $mt_gaz 'SCAL' 'CENTRE' Wga_gaz0 ;
 rvi . 'Wga_tfe' = 'KCHT' $mt_tfe 'SCAL' 'CENTRE' Wga_tfe0 ;
finsi ;
si egaliq ;
 rvi . 'Wgahe' = 'KCHT' $mthe 'SCAL' 'CENTRE' Wgahe0 ;
 rvi . 'Wgahe_p' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' Wgahe_p0 ;
finsi;
 
***Copy of the fields to do the convergence check :
thep = 'COPIER' (rvi . 'THE') ;
tcbp = 'COPIER' (rvi . 'TCB') ;
 
 
***Initialisation des proprietes physiques de l'eau :
rhe muhe lhe cphe = PRPEAU prt ro_o mu_o cp_o lamb_o theek ;
 
* Fluide chaud :
rvi . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhe ;
rvi . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhe cphe) ;
uhe = ELNO $mthe (KOPS qhe / rhe) ;
uhe = NOMC 'UZ' uhe ;
rvi . 'UHE' = 'KCHT' $mthe 'VECT' 'SOMMET' uhe ;
 
* Fluide froid :
rvi . 'LHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' lhe ;
rvi . 'RCPHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' ('*' rhe cphe) ;
uhe_p = ELNO $mthe_p (KOPS qhe_p / rhe) ;
uhe_p = NOMC 'UZ' (-1 * uhe_p) ;
rvi . 'UHE_p' = 'KCHT' $mthe_p 'VECT' 'SOMMET' uhe_p ;
 
 
fb1 = (exco uz (rvi . 'UHE')) ;
fbp1 = (exco uz (rvi . 'UHE_p')) ;
mess 'Initialisation : ';
mess ' -------------------------------' ;
mess 'rvi . uhe   : ' (mini fb1) (maxi fb1) ;
mess 'rvi . uhe_p : ' (mini fbp1) (maxi fbp1) ;
mess ' -------------------------------' ;
mess ;
 
 
***********FIN DE LA PREPARATION DU CALCUL MONOPHASIQUE*****************
************************************************************************
 
 
* BOUCLE DE CALCUL MONOPHASIQUE
************************************************************************
 
* Commentaires :
*  On affecte cette boucle d'un check de la convergence. Pour cela, on
*  cree un bouleen 'maxiter' qui permet d'arreter la boucle une fois que
*  la preision desiree sur la convergence est atteinte.
*  Le check de convergence permet aussi de prevenir l'utilisateur si le
*  systeme ne converge pas.
 
maxiter = vrai ;
 
'REPETER' bouc nb_bouc ;
MESS' 1) Calcul de la temperature moyenne dans l eau' ;
   tmoyhe = 'NOEL' $mthe (rvi . 'THE') ;
      'MESSAGE' 'Maxi Temperature moyenne fluide chaud (degres C) : '
             ('-' ('MAXIMUM' tmoyhe) CK) ;
   tmoyhe_p = 'NOEL' $mthe_p (rvi . 'THE') ;
      'MESSAGE' 'Maxi Temperature moyenne fluide froid (degres C) : '
             ('-' ('MAXIMUM' tmoyhe_p) CK) ;
 t_eau = ( REDU (rvi . 'THE') mthe ) ;
 t_eau_p = ( REDU (rvi . 'THE') mthe_p ) ;
tcavgaz = ( REDU (rvi . 'TCB') lcav31 ) ;
chpobor 'Tre fluide chaud en K' t_eau ;
chpobor 'Tre fluide froid en K' t_eau_p ;
chpobor 'Cavite gaz en K      ' tcavgaz ;
 
MESS '*  2) Transfert sur l interface';
   tmtif  = '*' mcf tmoyhe ; comm crayon fluide chaud ;
   tmtif1  = '*' mcf_q tmoyhe ; comm panier fluide chaud ;
   tmtif2  = '*' mcf_p tmoyhe_p ; comm panier fluide froid ;
   tmtif3  = '*' mcf_t tmoyhe_p ; comm t. force interne fluide froid ;
 
   rvi . 'T0' = 'KCHT' $mtif 'SCAL' 'CENTRE' tmtif ;
   rvi . 'T01' = 'KCHT' $mtin_p 'SCAL' 'CENTRE' tmtif1 ;
   rvi . 'T02' = 'KCHT' $mtif_pe 'SCAL' 'CENTRE' tmtif2 ;
   rvi . 'T03' = 'KCHT' $mtin_t 'SCAL' 'CENTRE' tmtif3;
 
MESS ' 3) Calcul du coefficient d echange a l interface   ' ;
*   3) Calcul du coefficient d'echange a l'interface
*      et des proprietes physiques du crayon chauffant
*      a partir des temperatures a l'instant precedent.
* Temperature de melange
   tm  = rvi . 'T0' ;
   tm1 = rvi . 'T01' ;
   tm2 = rvi . 'T02' ;
   tm3 = rvi . 'T03' ;
* Temperature a la paroi
   tp  = 'NOEL' $mtif   (rvi . 'TCB') ;
   tp1 = 'NOEL' $mtin_p (rvi . 'TCB') ;
   tp2 = 'NOEL' $mtif_pe (rvi . 'TCB') ;
   tp3 = 'NOEL' $mtin_t (rvi . 'TCB') ;
 
   t0 = '*' ('+' tm tp) 0.5D0 ;
   t01 = '*' ('+' tm1 tp1) 0.5D0 ;
   t02 = '*' ('+' tm2 tp2) 0.5D0 ;
   t03 = '*' ('+' tm3 tp3) 0.5D0 ;
 
   rho0 mu0 l0 cp0 = PRPEAU prt ro_o mu_o cp_o lamb_o t0 ;
   rho1 mu1 l1 cp1 = PRPEAU prt ro_o mu_o cp_o lamb_o t01 ;
   rho2 mu2 l2 cp2 = PRPEAU prt ro_o mu_o cp_o lamb_o t02 ;
   rho3 mu3 l3 cp3 = PRPEAU prt ro_o mu_o cp_o lamb_o t03 ;
   pran0 = '*' ('*' mu0 cp0) ('INVERSE' l0) ;
   pran1 = '*' ('*' mu1 cp1) ('INVERSE' l1) ;
   pran2 = '*' ('*' mu2 cp2) ('INVERSE' l2) ;
   pran3 = '*' ('*' mu3 cp3) ('INVERSE' l3) ;
 
* on prends la valeur absolue du debit s'il est < 0
  reyn0 = ( ( abs qhe_ga ) * diamh ) * ( INVE mu0 ) ;
  reyn1 = ( ( abs qhe_pa ) * diamh_p ) * ( INVE mu1 ) ;
  reyn2 = ( ( abs qhe_pe ) * diamh_pe ) * ( INVE mu2 ) ;
  reyn3 = ( ( abs qhe_t ) * diamh_t ) * ( INVE mu3 ) ;
 
* Correlation de Colburn
   nuss0 = A_colb '*' ('**' reyn0 B_colb)
                 '*' ('**' pran0 C_colb) ;
   nuss1 = A_colb '*' ('**' reyn1 B_colb)
                 '*' ('**' pran1 C_colb) ;
   nuss2 = A_colb '*' ('**' reyn2 B_colb)
                 '*' ('**' pran2 C_colb) ;
   nuss3 = A_colb '*' ('**' reyn3 B_colb)
                 '*' ('**' pran3 C_colb) ;
 
   h0 = nuss0 * l0 / diamh ;
   h1 = nuss1 * l1 / diamh_p ;
   h2 = nuss2 * l2 / diamh_pe ;
   h3 = nuss3 * l3 / diamh_t ;
 
****Test des coefficients d'echange h:
***h0 = h0 * 0.1;
***h1 = h1 * 50. ;
***h2 = h2 * 50. ;
***h3 = h3 * 100. ;
 
   rvi . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' h0 ;
   rvi . 'H1' =  'KCHT' $mtin_p 'SCAL' 'CENTRE' h1 ;
   rvi . 'H2' =  'KCHT' $mtif_pe 'SCAL' 'CENTRE' h2 ;
   rvi . 'H3' =  'KCHT' $mtin_t 'SCAL' 'CENTRE' h3 ;
 
* Proprietes physiques de la rallonge inferieure
   tcbckAr = 'NOEL' $mtcbAr (rvi . 'TCB') ;
   ju64 lcb = PRPZIRC tcbckAr ;
   rvi . 'LCBAr' = 'KCHT' $mtcbAr 'SCAL' 'CENTRE' lcb ;
 
* Proprietes physiques de la rallonge superieure
   tcbckEr = 'NOEL' $mtcbEr (rvi . 'TCB') ;
   ju64 lcb = PRPZIRC tcbckEr ;
   rvi . 'LCBEr' = 'KCHT' $mtcbEr 'SCAL' 'CENTRE' lcb ;
 
* Proprietes physiques du crayon chauffant
   tcbck = 'NOEL' $mtcb (rvi . 'TCB') ;
   ju64 lcb = PRPZIRC tcbck ;
   rvi . 'LCB' = 'KCHT' $mtcb 'SCAL' 'CENTRE' lcb ;
 
* Proprietes physiques du panier
   tcbck_p = 'NOEL' $mtcb_p (rvi . 'TCB') ;
   ju64 lcb_p = PRPZIRC tcbck_p ;
   rvi . 'LCB_p' = 'KCHT' $mtcb_p 'SCAL' 'CENTRE' lcb_p ;
 
* Proprietes physiques de la lame d'eau
   tck_lam = 'NOEL' $mt_lam (rvi . 'TCB') ;
   ju1 ju2 lcb_lam ju4 = PRPEAU prt ro_o mu_o cp_o lamb_o tck_lam ;
 
   rvi . 'LCB_lam' = 'KCHT' $mt_lam 'SCAL' 'CENTRE' lcb_lam ;
 
* Proprietes physiques du porte echantillon
   tck_pe = 'NOEL' $mt_pe (rvi . 'TCB') ;
   ju64 lcb_pe = PRPZIRC tck_pe ;
   rvi . 'LCB_pe' = 'KCHT' $mt_pe 'SCAL' 'CENTRE' lcb_pe ;
 
* Proprietes physiques du tube de force interne
   tcbck_i = 'NOEL' $mt_tfi (rvi . 'TCB') ;
   ju64 lcb_i = PRPZIRC tcbck_i ;
   rvi . 'LCB_tfi' = 'KCHT' $mt_tfi 'SCAL' 'CENTRE' lcb_i ;
 
* Proprietes physiques du tube de force externe
   tcbck_e = 'NOEL' $mt_tfe (rvi . 'TCB') ;
   ju64 lcb_e = PRPZIRC tcbck_e ;
   rvi . 'LCB_tfe' = 'KCHT' $mt_tfe 'SCAL' 'CENTRE' lcb_e ;
 
* Proprietes physiques du gaz entre tubes de force int. & ext.
   tcbck_g = 'NOEL' $mt_gaz (rvi . 'TCB') ;
   rhoju lcb_g  = PRPN2 tcbck_g ; comm diazote ;
   rvi . 'LCB_gaz' = 'KCHT' $mt_gaz 'SCAL' 'CENTRE' lcb_g ;
 
MESS
'3bis) Activation des sources internes dans les 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.
 
TKcbAr = NOEL $mtcbAr (rvi . 'TCB') ;
TKcbEr = NOEL $mtcbEr (rvi . 'TCB') ;
TKcb = NOEL $mtcb (rvi . 'TCB') ;
TKcb_p = NOEL $mtcb_p (rvi . 'TCB') ;
TK_lam = NOEL $mt_lam (rvi . 'TCB') ;
TK_pe = NOEL $mt_pe (rvi . 'TCB') ;
TK_tfi = NOEL $mt_tfi (rvi . 'TCB') ;
TK_gaz = NOEL $mt_gaz (rvi . 'TCB') ;
TK_tfe = NOEL $mt_tfe (rvi . 'TCB') ;
 
rhocbAr lju = PRPZIRC TKcbAr ;
rhocbEr lju = PRPZIRC TKcbEr ;
rhocb lju = PRPZIRC TKcb ;
rhocb_p lju = PRPZIRC TKcb_p ;
rho_lam lju = PRPZIRC TK_lam ;
rho_pe lju = PRPZIRC TK_pe ;
rho_tfi lju = PRPZIRC TK_tfi ;
rho_gaz lju = PRPN2 TK_gaz ;
rho_tfe lju = PRPZIRC TK_tfe ;
 
si egasol ;
 rvi . 'WgacbAr' = KOPS rhocbAr * EgacbAr ;
 rvi . 'WgacbEr' = KOPS rhocbEr * EgacbEr ;
 rvi . 'Wgacb' = KOPS rhocb * Egacb ;
 rvi . 'Wgacb_p' = KOPS rhocb_p * Egacb_p ;
 rvi . 'Wga_lam' = KOPS rho_lam * Ega_lam ;
 rvi . 'Wga_pe' = KOPS rho_pe * Ega_pe ;
 rvi . 'Wga_tfi' = KOPS rho_tfi * Ega_tfi ;
 rvi . 'Wga_gaz' = KOPS rho_gaz * Ega_gaz ;
 rvi . 'Wga_tfe' = KOPS rho_tfe * Ega_tfe ;
finsi ;
 
MESS '3)ter Condition de convection naturelle a la piscine' ;
*   3)ter Condition de convection naturelle a la piscine
Tef_t = NOEL $mtef_t (rvi . 'TCB') ;
SI varcna ;
mess 'Convection naturelle avec la piscine' ;
 h_pisci = HNAT Tef_t tpiscik (zM - z1) ;
SINON ;
mess 'Echange constant avec la piscine : ' hcncst ' W/m2.K' ;
 h_pisci = hcncst ;
FINSI ;
rvi . 'Hn4' = KCHT $mtef_t 'SCAL' 'CENTRE' h_pisci ;
mess ( mini (rvi . 'Hn4')) ( maxi (rvi . 'Hn4')) ;
SI rayo  ;
 
MESS '3 qua) Rayonnement face a face' ;
*   3 qua) Rayonnement face a face'
t2i_ = ( REDU (rvi . 'TCB') dcav_2 ) ;
t2i_pro = KPRO t2i_ ( INVE rel_12 ) ;
rvcb . inco . 'T2I' = KCHT $mtif_t SCAL SOMMET t2i_pro ; comm rayoju ;
 
t1s_ = ( REDU (rvi . 'TCB') dcav_1 ) ;
t1s_pro = KPRO t1s_  rel_12  ;
rvcb . inco . 'T1S' = KCHT $mten_t SCAL SOMMET t1s_pro ; comm rayoju ;
 
2t1 = t1s_*t1s_  ;
2t2 = t2i_pro * t2i_pro ;
h_1 = stefan * alpha1 * (2t1 + 2t2) * (t1s_ + t2i_pro) ;
 
h_1c = NOEL $mtif_t h_1 ; comm rayoju ;
 
h_2  = KPRO h_1 rel_12 ;
h_2c = NOEL $mten_t h_2 ; comm rayoju ;
 
rvcb . inco . 'HR1' = KCHT $mtif_t SCAL CENTRE h_1c ; comm rayoju ;
rvcb . inco . 'HR2' = KCHT $mten_t SCAL CENTRE h_2c ; comm rayoju ;
FINSI ; comm rayonnement face a face ;
 
 
MESS ' 4) Resolution de la thermique dans le crayon chauffant ' ;
*   4) Resolution de la thermique dans le crayon chauffant
   EXEC rvcb ;
 
 
MESS  '  5) Calcul du flux de chaleur aux interfaces Solides/Fluides' ;
*   5) Calcul du flux de chaleur aux interfaces Solides/Fluides
 
***Gaine/Fluide chaud
   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 du crayon (W) : '
            ('MAXIMUM' ('RESULT' f0)) ;
 
***Panier/Fluide chaud
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_p 'OPER' 'MDIA' 'H1' 'INCO' 'TCB' 'F01' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_p 'OPER' 'MDIA' 'H1' 'INCO' 'T01'  'F01' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f01 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T01'  (rvi . 'T01')))) ;
        'MESSAGE' 'Flux integre sur la paroi chaude du panier (W) : '
             ('MAXIMUM' ('RESULT' f01)) ;
 
***Panier/Fluide froid
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif_pe 'OPER' 'MDIA' 'H2' 'INCO' 'TCB' 'F02' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif_pe 'OPER' 'MDIA' 'H2' 'INCO' 'T02'  'F02' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f02 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T02'  (rvi . 'T02')))) ;
      'MESSAGE' 'Flux integre sur la paroi froide du PE (W) : '
             ('MAXIMUM' ('RESULT' f02)) ;
 
***Tube de force/Fluide froid
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_t 'OPER' 'MDIA' 'H3' 'INCO' 'TCB' 'F03' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_t 'OPER' 'MDIA' 'H3' 'INCO' 'T03'  'F03' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f03 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T03'  (rvi . 'T03')))) ;
      'MESSAGE'
  'Flux integre sur la paroi interne du tube de force interne (W) : '
             ('MAXIMUM' ('RESULT' f03)) ;
 
***Tube de force /Piscine
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtef_t 'OPER' 'MDIA' 'Hn4' 'INCO' 'TCB' 'Fn4' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtef_t 'OPER' 'MDIA' 'Hn4' 'INCO' 'Tp4'  'Fn4' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   fn4 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'Tp4'  (rvi . 'Tp4')))) ;
      'MESSAGE'
  'Flux integre sur la paroi externe du tube de force externe (W) : '
             ('MAXIMUM' ('RESULT' fn4)) ;
 
 
*   6) Transfert en une puissance volumique dans l'eau
   rvi . 'F0' = f0 ; comm interface crayon et fluide chaud ;
   phe = '*' mfc f0 ;
   rvi . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe ;
 
   rvi . 'F02' = f02 ; comm interface panier et fluide froid ;
   phe2 = '*' mfc_p f02 ;
   rvi . 'PHE2' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' phe2 ;
 
   rvi . 'F01' = f01 ; comm interface panier et fluide chaud ;
   phe1 = '*' mfc_q f01 ;
   rvi . 'PHE1' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe1 ;
 
   rvi . 'F03' = f03 ; comm interface t. force et fluide froid ;
   phe3 = '*' mfc_t f03 ;
   rvi . 'PHE3' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' phe3 ;
 
mess 'Puissance volumique fournie a l eau : ' ;
mess 'interface crayon et fluide chaud   : '(MINI phe) (MAXi phe) ;
mess 'interface panier et fluide chaud   : '(MINI phe1) (MAXi phe1) ;
mess 'interface panier et fluide froid   : '(MINI phe2) (MAXi phe2) ;
mess 'interface t. force et fluide froid : '(MINI phe3) (MAXi phe3) ;
 
MESS '6bis) Activation des sources internes dans l eau (en W/m3)' ;
*   6bis) Activation des sources internes dans de l'eau (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 dans.
 
rhoeauc ju1 ju2 ju3 = PRPEAU prt ro_o mu_o cp_o lamb_o tmoyhe ;
rhoeauf ju1 ju2 ju3 = PRPEAU prt ro_o mu_o cp_o lamb_o tmoyhe_p ;
 
si egaliq ;
 rvi . 'Wgahe' = rhoeauc * Egahe ;
 rvi . 'Wgahe_p' = rhoeauf * Egahe_p ;
finsi ;
 
MESS ' 7) Calcul des proprietes physiques de l eau ' ;
*   7) Calcul des proprietes physiques de l'eau
*      a partir des variables a l'instant precedent.
*
* Fluide chaud :
dmthe = DOMA $mthe 'MAILLAGE' ;
Thes = REDU dmthe (rvi . 'THE') ;
Thec = 'NOEL' $mthe (rvi . 'THE') ;
rhes muhes lhes cphes = PRPEAU prt ro_o mu_o cp_o lamb_o Thes ;
rhec muhec lhec cphec = PRPEAU prt ro_o mu_o cp_o lamb_o Thec ;
uhes = ELNO $mthe (KOPS qhe / rhec) ;
uhes = NOMC 'UZ' uhes ;
rvi . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhec cphec) ;
rvi . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhec ;
rvi . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhes ;
 
* Fluide froid :
dmthe_p = DOMA $mthe_p 'MAILLAGE' ;
Thes = REDU dmthe_p (rvi . 'THE') ;
Thec = 'NOEL' $mthe_p (rvi . 'THE') ;
rhes muhes lhes cphes = PRPEAU prt ro_o mu_o cp_o lamb_o Thes ;
rhec muhec lhec cphec = PRPEAU prt ro_o mu_o cp_o lamb_o Thec ;
uhes = ELNO $mthe_p (KOPS qhe_p / rhec) ;
uhes = NOMC 'UZ' uhes ;
uhes = NOMC 'UZ' (-1 * uhes) ;
rvi . 'RCPHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' ('*' rhec cphec) ;
rvi . 'LHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' lhec ;
rvi . 'UHE_p' =  'KCHT' $mthe_p 'VECT' 'SOMMET' uhes ;
 
 
fb1 = (exco uz (rvi . 'UHE')) ;
fbp1 = (exco uz (rvi . 'UHE_p')) ;
mess ;
mess ' -------------------------------' ;
mess 'rvi . uhe   : ' (mini fb1) (maxi fb1) ;
mess 'rvi . uhe_p : ' (mini fbp1) (maxi fbp1) ;
mess ' -------------------------------' ;
mess ;
 
 
MESS '  8) Calcul de la temperature dans l eau' ;
*   8) Calcul de la temperature dans l'eau
   EXEC rvhe ;
 
***Check de Convergence :
   then = 'COPIER' (rvi . 'THE') ;
   tcbn = 'COPIER' (rvi . 'TCB') ;
   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 ;
****
 
'SI' ( 'ET' ('EGA' complet vrai) ('EGA' maxiter vrai) );
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'The computation has not converged') ;
   'ERREUR' 'Check calisa1.procedur !' ;
'FINSI' ;
 
***********FIN DE LA BOUCLE DE CALCUL MONOPHASIQUE**********************
************************************************************************
 
tcpu9 ;
 
*****------------PASSAGE EN DIPHASIQUE SI NECESSAIRE------------********
************************************************************************
 
si ( NON complet ) ;
  ebulli1 = FAUX ;
sinon ;
 ebulli1 c_mono1 c_diph1 c_mono1N c_diph1N = EBULL ;
finsi ;
 
SI (NON ebulli1) ;
   MESS 'On reste en monophasique , pas d ebullition nuclee' ;
   MESS 'On peut passer au post traitement' ;
   opti echo 0 ;
 
SINON ;
   MESS 'ebullition locale sur la gaine, c un probleme diphasique' ;
   MESS 'il faut passer au jeu de donnees diphasique modifiee' ;
 
* PREPARATION DU CALCUL DIPHASIQUE de "essai11"
************************************************************************
 
*      IV.3. Initialisation de toutes les variables avant le calcul
rvi . 'TCB' = 'KCHT'
    ( $mtcbAr et $mtcbEr et $mtcb et $mtcb_p et $mt_lam et $mt_pe et
      $mt_tfi et $mt_tfe et $mt_gaz ) 'SCAL' 'SOMMET' theek ;
 
ju64 lcb = PRPZIRC theek ;
rhoju lcb_g = PRPN2 theek ; comm diazote ;
**Inox a rajouter
 
rvi . 'LCBAr' = 'KCHT' $mtcbAr 'SCAL' 'CENTRE' lcb ;
rvi . 'LCBEr' = 'KCHT' $mtcbEr 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB' = 'KCHT' $mtcb 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_p' = 'KCHT' $mtcb_p 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_lam' = 'KCHT' $mt_lam 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_pe' = 'KCHT' $mt_pe 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_tfi' = 'KCHT' $mt_tfi 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_tfe' = 'KCHT' $mt_tfe 'SCAL' 'CENTRE' lcb ;
rvi . 'LCB_gaz' = 'KCHT' $mt_gaz 'SCAL' 'CENTRE' lcb_g ;
 
rvi . 'PCB' = 'KCHT' $mtin 'SCAL' 'CENTRE' pmsurf ;
rvi . 'T0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T01' = 'KCHT'  $mtin_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T02' = 'KCHT'  $mtif_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T03' = 'KCHT'  $mtin_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F01' = 'KCHT'  $mtin_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F02' = 'KCHT'  $mtif_pe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'F03' = 'KCHT'  $mtin_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H'  = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H1'  = 'KCHT'  $mtin_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H2'  = 'KCHT'  $mtif_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'H3'  = 'KCHT'  $mtin_t 'SCAL' 'CENTRE' 0.D0 ;
SI rayo ;
rvi . 'T1S' = 'KCHT'  $mten_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'T2I' = 'KCHT'  $mtif_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'HR2' = 'KCHT'  $mten_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'HR1' = 'KCHT'  $mtif_t 'SCAL' 'CENTRE' 0.D0 ;
FINSI ; comm rayonnement face a face ;
rvi . 'THE' = 'KCHT' ( $mthe ET $mthe_p ) 'SCAL' 'SOMMET' theek ;
rvi . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PHE1' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PHE2' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'PHE3' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' 0.D0 ;
 
*Essai27 : changement de la CL a la piscine :
*------------------------------------------------
rvi . 'Hn4'  = 'KCHT' $mtef_t 'SCAL' 'CENTRE' 0.D0 ;
rvi . 'Tp4' = 'KCHT' $mtef_t 'SCAL' 'CENTRE' tpiscik ;
rvi . 'Fn4' = 'KCHT' $mtef_t 'SCAL' 'CENTRE' 0.D0 ;
*------------------------------------------------
 
***Initialisation des sources internes :
si egasol ;
 rvi . 'WgacbAr' = 'KCHT' $mtcbAr 'SCAL' 'CENTRE' WgacbAr0 ;
 rvi . 'WgacbEr' = 'KCHT' $mtcbEr 'SCAL' 'CENTRE' WgacbEr0 ;
 rvi . 'Wgacb' = 'KCHT' $mtcb 'SCAL' 'CENTRE' Wgacb0 ;
 rvi . 'Wgacb_p' = 'KCHT' $mtcb_p 'SCAL' 'CENTRE' Wgacb_p0 ;
 rvi . 'Wga_lam' = 'KCHT' $mt_lam 'SCAL' 'CENTRE' Wga_lam0 ;
 rvi . 'Wga_pe' = 'KCHT' $mt_pe 'SCAL' 'CENTRE' Wga_pe0 ;
 rvi . 'Wga_tfi' = 'KCHT' $mt_tfi 'SCAL' 'CENTRE' Wga_tfi0 ;
 rvi . 'Wga_gaz' = 'KCHT' $mt_gaz 'SCAL' 'CENTRE' Wga_gaz0 ;
 rvi . 'Wga_tfe' = 'KCHT' $mt_tfe 'SCAL' 'CENTRE' Wga_tfe0 ;
finsi ;
si egaliq ;
 rvi . 'Wgahe' = 'KCHT' $mthe 'SCAL' 'CENTRE' Wgahe0 ;
 rvi . 'Wgahe_p' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' Wgahe_p0 ;
finsi ;
 
Tpsat = KCHT $mtif 'SCAL' 'CENTRE' Tsat ;
f0 = KCHT $mtif 'SCAL' 'CENTRE' 0. ;
 
***Copy of the fields to do the convergence check :
thep = 'COPIER' (rvi . 'THE') ;
tcbp = 'COPIER' (rvi . 'TCB') ;
 
***Initialisation des proprietes physiques de l'eau :
rhe muhe lhe cphe = PRPEAU prt ro_o mu_o cp_o lamb_o theek ;
 
* Fluide chaud :
rvi . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhe ;
rvi . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhe cphe) ;
uhe = ELNO $mthe (KOPS qhe / rhe) ;
uhe = NOMC 'UZ' uhe ;
rvi . 'UHE' = 'KCHT' $mthe 'VECT' 'SOMMET' uhe ;
 
* Fluide froid :
rvi . 'LHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' lhe ;
rvi . 'RCPHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' ('*' rhe cphe) ;
uhe_p = ELNO $mthe_p (KOPS qhe_p / rhe) ;
uhe_p = NOMC 'UZ' (-1 * uhe_p) ;
rvi . 'UHE_p' = 'KCHT' $mthe_p 'VECT' 'SOMMET' uhe_p ;
 
 
fb1 = (exco uz (rvi . 'UHE')) ;
fbp1 = (exco uz (rvi . 'UHE_p')) ;
mess 'Initialisation : ';
mess ' -------------------------------' ;
mess 'rvi . uhe   : ' (mini fb1) (maxi fb1) ;
mess 'rvi . uhe_p : ' (mini fbp1) (maxi fbp1) ;
mess ' -------------------------------' ;
mess ;
 
 
*******FIN DE LA PREPARATION DU CALCUL DIPHASIQUE de "essai11"**********
************************************************************************
 
 
* BOUCLE DE CALCUL DIPHASIQUE de "essai11"
************************************************************************
 
* Commentaires :
*  On affecte cette boucle d'un check de la convergence. Pour cela, on
*  cree un bouleen 'maxiter' qui permet d'arreter la boucle une fois que
*  la preision desiree sur la convergence est atteinte.
*  Le check de convergence permet aussi de prevenir l'utilisateur si le
*  systeme ne converge pas.
 
maxiter = vrai ;
 
'REPETER' bouc nb_bouc ;
MESS' 1) Calcul de la temperature moyenne dans l eau' ;
   tmoyhe = 'NOEL' $mthe (rvi . 'THE') ;
      'MESSAGE' 'Maxi Temperature moyenne fluide chaud (degres C) : '
             ('-' ('MAXIMUM' tmoyhe) CK) ;
   tmoyhe_p = 'NOEL' $mthe_p (rvi . 'THE') ;
      'MESSAGE' 'Maxi Temperature moyenne fluide froid (degres C) : '
             ('-' ('MAXIMUM' tmoyhe_p) CK) ;
 t_eau = ( REDU (rvi . 'THE') mthe ) ;
 t_eau_p = ( REDU (rvi . 'THE') mthe_p ) ;
 tcavgaz = ( REDU (rvi . 'TCB') lcav31 ) ;
chpobor 'Tre fluide chaud en K' t_eau ;
chpobor 'Tre fluide froid en K' t_eau_p ;
chpobor 'Cavite gaz en K      ' tcavgaz ;
 
MESS '*  2) Transfert sur l interface';
   tmtif  = '*' mcf tmoyhe ; comm crayon fluide chaud ;
   tmtif1  = '*' mcf_q tmoyhe ; comm panier fluide chaud ;
   tmtif2  = '*' mcf_p tmoyhe_p ; comm panier fluide froid ;
   tmtif3  = '*' mcf_t tmoyhe_p ; comm t. force interne fluide froid ;
 
****RAJOUT ESSAI14**** :
***Diphasique : determination de la temperature de saturation a la paroi
* Remarque :
*   Dans la formule de Jens et Lottes, il faut faire attention a prendre
*   la pression en bars en le flux en W/cm2
 
* Determination de la temperature de paroi a saturation :
cf0 = f0 MASQ 'SUPERIEUR' 0. ;
f0 = KOPS f0 * cf0 ;
Volif = DOMA $mtif 'VOLUME' ;
 phi0 = KOPS f0 / Volif;
 phi0 = 1.e-4 * phi0 ;
 
expo1 = EXP ( -1. * PsatB / 62.) ;
 
DTsat = 7.91 * expo1 *(phi0 ** 0.25) ;
Tpsat = Tsat + DTsat ;
****
 
T0_md = (c_mono1 * tmtif) + (c_diph1 * Tpsat) ; comm modif essai14 ;
 
   rvi . 'T0' = 'KCHT' $mtif 'SCAL' 'CENTRE' T0_md ; comm mod essai14;
   rvi . 'T01' = 'KCHT' $mtin_p 'SCAL' 'CENTRE' tmtif1 ;
   rvi . 'T02' = 'KCHT' $mtif_pe 'SCAL' 'CENTRE' tmtif2;
   rvi . 'T03' = 'KCHT' $mtin_t 'SCAL' 'CENTRE' tmtif3;
 
MESS ' 3) Calcul du coefficient d echange a l interface   ' ;
*   3) Calcul du coefficient d'echange a l'interface
*      et des proprietes physiques du crayon chauffant
*      a partir des temperatures a l'instant precedent.
* Temperature de melange
   tm  = rvi . 'T0' ;
   tm1 = rvi . 'T01' ;
   tm2 = rvi . 'T02' ;
   tm3 = rvi . 'T03' ;
* Temperature a la paroi
   tp  = 'NOEL' $mtif   (rvi . 'TCB') ;
   tp1 = 'NOEL' $mtin_p (rvi . 'TCB') ;
   tp2 = 'NOEL' $mtif_pe (rvi . 'TCB') ;
   tp3 = 'NOEL' $mtin_t (rvi . 'TCB') ;
 
   t0 = '*' ('+' tm tp) 0.5D0 ;
   t01 = '*' ('+' tm1 tp1) 0.5D0 ;
   t02 = '*' ('+' tm2 tp2) 0.5D0 ;
   t03 = '*' ('+' tm3 tp3) 0.5D0 ;
 
   rho0 mu0 l0 cp0 = PRPEAU prt ro_o mu_o cp_o lamb_o t0 ;
   rho1 mu1 l1 cp1 = PRPEAU prt ro_o mu_o cp_o lamb_o t01 ;
   rho2 mu2 l2 cp2 = PRPEAU prt ro_o mu_o cp_o lamb_o t02 ;
   rho3 mu3 l3 cp3 = PRPEAU prt ro_o mu_o cp_o lamb_o t03 ;
   pran0 = '*' ('*' mu0 cp0) ('INVERSE' l0) ;
   pran1 = '*' ('*' mu1 cp1) ('INVERSE' l1) ;
   pran2 = '*' ('*' mu2 cp2) ('INVERSE' l2) ;
   pran3 = '*' ('*' mu3 cp3) ('INVERSE' l3) ;
 
* on prends la valeur absolue du debit s'il est < 0
  reyn0 = ( ( abs qhe_ga ) * diamh ) * ( INVE mu0 ) ;
  reyn1 = ( ( abs qhe_pa ) * diamh_p ) * ( INVE mu1 ) ;
  reyn2 = ( ( abs qhe_pe ) * diamh_pe ) * ( INVE mu2 ) ;
  reyn3 = ( ( abs qhe_t ) * diamh_t ) * ( INVE mu3 ) ;
 
* Correlation de Colburn
   nuss0 = A_colb '*' ('**' reyn0 B_colb)
                 '*' ('**' pran0 C_colb) ;
   nuss1 = A_colb '*' ('**' reyn1 B_colb)
                 '*' ('**' pran1 C_colb) ;
   nuss2 = A_colb '*' ('**' reyn2 B_colb)
                 '*' ('**' pran2 C_colb) ;
   nuss3 = A_colb '*' ('**' reyn3 B_colb)
                 '*' ('**' pran3 C_colb) ;
 
   h0_m = nuss0 * l0 / diamh ;
   h0_d = 1.e12 ;
   h0 = (c_mono1 * h0_m) + (c_diph1 * h0_d) ; comm modif essai14 ;
 
   h1 = nuss1 * l1 / diamh_p ;
   h2 = nuss2 * l2 / diamh_pe ;
   h3 = nuss3 * l3 / diamh_t ;
 
****Test des coefficients d'echange h:
***h0 = h0 * 0.1;
***h1 = h1 * 50. ;
***h2 = h2 * 50. ;
***h3 = h3 * 1.e-8 ;
 
   rvi . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' h0 ;
   rvi . 'H1' =  'KCHT' $mtin_p 'SCAL' 'CENTRE' h1 ;
   rvi . 'H2' =  'KCHT' $mtif_pe 'SCAL' 'CENTRE' h2 ;
   rvi . 'H3' =  'KCHT' $mtin_t 'SCAL' 'CENTRE' h3 ;
 
* Proprietes physiques de la rallonge inferieure
   tcbckAr = 'NOEL' $mtcbAr (rvi . 'TCB') ;
   ju64 lcb = PRPZIRC tcbckAr ;
   rvi . 'LCBAr' = 'KCHT' $mtcbAr 'SCAL' 'CENTRE' lcb ;
 
* Proprietes physiques de la rallonge superieure
   tcbckEr = 'NOEL' $mtcbEr (rvi . 'TCB') ;
   ju64 lcb = PRPZIRC tcbckEr ;
   rvi . 'LCBEr' = 'KCHT' $mtcbEr 'SCAL' 'CENTRE' lcb ;
 
* Proprietes physiques du crayon chauffant
   tcbck = 'NOEL' $mtcb (rvi . 'TCB') ;
   ju64 lcb = PRPZIRC tcbck ;
   rvi . 'LCB' = 'KCHT' $mtcb 'SCAL' 'CENTRE' lcb ;
 
* Proprietes physiques du panier
   tcbck_p = 'NOEL' $mtcb_p (rvi . 'TCB') ;
   ju64 lcb_p = PRPZIRC tcbck_p ;
   rvi . 'LCB_p' = 'KCHT' $mtcb_p 'SCAL' 'CENTRE' lcb_p ;
 
* Proprietes physiques de la lame d'eau
   tck_lam = 'NOEL' $mt_lam (rvi . 'TCB') ;
   ju1 ju2 lcb_lam ju4 = PRPEAU prt ro_o mu_o cp_o lamb_o tck_lam ;
   rvi . 'LCB_lam' = 'KCHT' $mt_lam 'SCAL' 'CENTRE' lcb_lam ;
 
* Proprietes physiques du porte echantillon
   tck_pe = 'NOEL' $mt_pe (rvi . 'TCB') ;
   ju64 lcb_pe = PRPZIRC tck_pe ;
   rvi . 'LCB_pe' = 'KCHT' $mt_pe 'SCAL' 'CENTRE' lcb_pe ;
 
* Proprietes physiques du tube de force interne
   tcbck_i = 'NOEL' $mt_tfi (rvi . 'TCB') ;
   ju64 lcb_i = PRPZIRC tcbck_i ;
   rvi . 'LCB_tfi' = 'KCHT' $mt_tfi 'SCAL' 'CENTRE' lcb_i ;
 
* Proprietes physiques du tube de force externe
   tcbck_e = 'NOEL' $mt_tfe (rvi . 'TCB') ;
   ju64 lcb_e = PRPZIRC tcbck_e ;
   rvi . 'LCB_tfe' = 'KCHT' $mt_tfe 'SCAL' 'CENTRE' lcb_e ;
 
* Proprietes physiques du gaz entre tubes de force int. & ext.
   tcbck_g = 'NOEL' $mt_gaz (rvi . 'TCB') ;
   rhoju lcb_g = PRPN2 tcbck_g ; comm diazote ;
   rvi . 'LCB_gaz' = 'KCHT' $mt_gaz 'SCAL' 'CENTRE' lcb_g ;
 
MESS
'3bis) Activation des sources internes dans les structures (en W/m3)' ;
*   3bis) Activation des sources internes dans les 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.
 
TKcbAr = NOEL $mtcbAr (rvi . 'TCB') ;
TKcbEr = NOEL $mtcbEr (rvi . 'TCB') ;
TKcb = NOEL $mtcb (rvi . 'TCB') ;
TKcb_p = NOEL $mtcb_p (rvi . 'TCB') ;
TK_lam = NOEL $mt_lam (rvi . 'TCB') ;
TK_pe = NOEL $mt_pe (rvi . 'TCB') ;
TK_tfi = NOEL $mt_tfi (rvi . 'TCB') ;
TK_gaz = NOEL $mt_gaz (rvi . 'TCB') ;
TK_tfe = NOEL $mt_tfe (rvi . 'TCB') ;
 
rhocbAr lju = PRPZIRC TKcbAr ;
rhocbEr lju = PRPZIRC TKcbEr ;
rhocb lju = PRPZIRC TKcb ;
rhocb_p lju = PRPZIRC TKcb_p ;
rho_lam lju = PRPZIRC TK_lam ;
rho_pe lju = PRPZIRC TK_pe ;
rho_tfi lju = PRPZIRC TK_tfi ;
rho_gaz lju = PRPN2 TK_gaz ;
rho_tfe lju = PRPZIRC TK_tfe ;
 
si egasol ;
 rvi . 'WgacbAr' = KOPS rhocbAr * EgacbAr ;
 rvi . 'WgacbEr' = KOPS rhocbEr * EgacbEr ;
 rvi . 'Wgacb' = KOPS rhocb * Egacb ;
 rvi . 'Wgacb_p' = KOPS rhocb_p * Egacb_p ;
 rvi . 'Wga_lam' = KOPS rho_lam * Ega_lam ;
 rvi . 'Wga_pe' = KOPS rho_pe * Ega_pe ;
 rvi . 'Wga_tfi' = KOPS rho_tfi * Ega_tfi ;
 rvi . 'Wga_gaz' = KOPS rho_gaz * Ega_gaz ;
 rvi . 'Wga_tfe' = KOPS rho_tfe * Ega_tfe ;
finsi ;
 
MESS '3)ter Condition de convection naturelle a la piscine' ;
*   3)ter Condition de convection naturelle a la piscine
Tef_t = NOEL $mtef_t (rvi . 'TCB') ;
 
SI varcna ;
mess 'Convection naturelle avec la piscine' ;
 h_pisci = HNAT Tef_t tpiscik (zM - z1) ;
SINON ;
mess 'Echange constant avec la piscine : ' hcncst ' W/m2.K' ;
 h_pisci = hcncst ;
FINSI ;
rvi . 'Hn4' = KCHT $mtef_t 'SCAL' 'CENTRE' h_pisci ;
mess ( mini (rvi . 'Hn4')) ( maxi (rvi . 'Hn4')) ;
 
bool1 = EGA &bouc 1 ;
SI rayo ;
 
MESS '3 qua) Rayonnement face a face' ;
* 3 qua) Rayonnement face a face'
t2i_ = ( REDU (rvi . 'TCB') dcav_2 ) ;
t2i_pro = KPRO t2i_ ( INVE rel_12 ) ;
rvcb . inco . 'T2I' = KCHT $mtif_t SCAL SOMMET t2i_pro ; comm rayoju ;
 
t1s_ = ( REDU (rvi . 'TCB') dcav_1 ) ;
t1s_pro = KPRO t1s_  rel_12  ;
rvcb . inco . 'T1S' = KCHT $mten_t SCAL SOMMET t1s_pro ; comm rayoju ;
 
2t1 = t1s_*t1s_  ;
2t2 = t2i_pro * t2i_pro ;
h_1 = stefan * alpha1 * (2t1 + 2t2) * (t1s_ + t2i_pro) ;
 
h_1c = NOEL $mtif_t h_1 ; comm rayoju ;
 
h_2  = KPRO h_1 rel_12 ;
h_2c = NOEL $mten_t h_2 ; comm rayoju ;
 
rvcb . inco . 'HR1' = KCHT $mtif_t SCAL CENTRE h_1c ; comm rayoju ;
rvcb . inco . 'HR2' = KCHT $mten_t SCAL CENTRE h_2c ; comm rayoju ;
FINSI ; comm rayonnement face a face ;
 
 
MESS ' 4) Resolution de la thermique dans le crayon chauffant ' ;
*   4) Resolution de la thermique dans le crayon chauffant
   EXEC rvcb ;
 
 
MESS
 '  5) Calcul du flux de chaleur a la surface du crayon chauffant' ;
*   5) Calcul du flux de chaleur a la surface du crayon chauffant
 
***Diphasique : creation d'un flux de chaleur conductif
     lambd1 = rvi . 'LCB' ;
          T = rvi . 'TCB' ;
Wif f0_d Flif = FLUXNS choiju lambd1 T $mtif $mtcb ;
 
SI (choiju EGA 1) ;
         f0_d = NOEL $mtif f0_d ;
FINSI;
 
***Monophasique : flux de chaleur de convection
   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_m = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                           ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                           ('NOMC' 'T0'  (rvi . 'T0')))) ;
 
 
f0 = (c_mono1 * f0_m) + (c_diph1 * f0_d) ; comm modif essai14 ;
 
     'MESSAGE' 'Flux integre sur la paroi du crayon (W) : '
            ('MAXIMUM' ('RESULT' f0)) ;
 
 
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_p 'OPER' 'MDIA' 'H1' 'INCO' 'TCB' 'F01' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_p 'OPER' 'MDIA' 'H1' 'INCO' 'T01'  'F01' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f01 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T01'  (rvi . 'T01')))) ;
        'MESSAGE' 'Flux integre sur la paroi chaude du panier (W) : '
             ('MAXIMUM' ('RESULT' f01)) ;
 
 
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif_pe 'OPER' 'MDIA' 'H2' 'INCO' 'TCB' 'F02' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif_pe 'OPER' 'MDIA' 'H2' 'INCO' 'T02'  'F02' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f02 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T02'  (rvi . 'T02')))) ;
      'MESSAGE' 'Flux integre sur la paroi froide du PE (W) : '
             ('MAXIMUM' ('RESULT' f02)) ;
 
 
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_t 'OPER' 'MDIA' 'H3' 'INCO' 'TCB' 'F03' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtin_t 'OPER' 'MDIA' 'H3' 'INCO' 'T03'  'F03' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f03 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T03'  (rvi . 'T03')))) ;
      'MESSAGE'
  'Flux integre sur la paroi interne du tube de force interne (W) : '
             ('MAXIMUM' ('RESULT' f03)) ;
 
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtef_t 'OPER' 'MDIA' 'Hn4' 'INCO' 'TCB' 'Fn4' ;
   rvmdia1 . 'INCO' = rvi ;
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtef_t 'OPER' 'MDIA' 'Hn4' 'INCO' 'Tp4'  'Fn4' ;
   rvmdia2 . 'INCO' = rvi ;
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   fn4 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCB' (rvi . 'TCB')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'Tp4'  (rvi . 'Tp4')))) ;
      'MESSAGE'
  'Flux integre sur la paroi externe du tube de force externe (W) : '
             ('MAXIMUM' ('RESULT' fn4)) ;
 
*   6) Transfert en une puissance volumique dans l'eau
   rvi . 'F0' = f0 ; comm interface crayon et fluide chaud ;
   phe = '*' mfc f0 ;
   rvi . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe ;
 
   rvi . 'F02' = f02 ; comm interface panier et fluide froid ;
   phe2 = '*' mfc_p f02 ;
   rvi . 'PHE2' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' phe2 ;
 
   rvi . 'F01' = f01 ; comm interface panier et fluide chaud ;
   phe1 = '*' mfc_q f01 ;
   rvi . 'PHE1' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe1 ;
 
   rvi . 'F03' = f03 ; comm interface t. force et fluide froid ;
   phe3 = '*' mfc_t f03 ;
   rvi . 'PHE3' = 'KCHT' $mthe_p 'SCAL' 'CENTRE' phe3 ;
 
mess 'Puissance volumique fournie a l eau : ' ;
mess 'interface crayon et fluide chaud   : '(MINI phe) (MAXi phe) ;
mess 'interface panier et fluide chaud   : '(MINI phe1) (MAXi phe1) ;
mess 'interface panier et fluide froid   : '(MINI phe2) (MAXi phe2) ;
mess 'interface t. force et fluide froid : '(MINI phe3) (MAXi phe3) ;
 
MESS
'6bis) Activation des sources internes dans l eau (en W/m3)' ;
*   6bis) Activation des sources internes de l'eau (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.
 
rhoeauc ju1 ju2 ju3 = PRPEAU prt ro_o mu_o cp_o lamb_o tmoyhe ;
rhoeauf ju1 ju2 ju3 = PRPEAU prt ro_o mu_o cp_o lamb_o tmoyhe_p ;
 
si egaliq ;
 rvi . 'Wgahe' = rhoeauc * Egahe ;
 rvi . 'Wgahe_p' = rhoeauf * Egahe_p ;
finsi ;
 
MESS ' 7) Calcul des proprietes physiques de l eau ' ;
*   7) Calcul des proprietes physiques de l'eau
*      a partir des variables a l'instant precedent.
*
* Fluide chaud :
dmthe = DOMA $mthe 'MAILLAGE' ;
Thes = REDU dmthe (rvi . 'THE') ;
Thec = 'NOEL' $mthe (rvi . 'THE') ;
rhes muhes lhes cphes = PRPEAU prt ro_o mu_o cp_o lamb_o Thes ;
rhec muhec lhec cphec = PRPEAU prt ro_o mu_o cp_o lamb_o Thec ;
uhes = ELNO $mthe (KOPS qhe / rhec) ;
uhes = NOMC 'UZ' uhes ;
rvi . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhec cphec) ;
rvi . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhec ;
rvi . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhes ;
 
 
* Fluide froid :
dmthe_p = DOMA $mthe_p 'MAILLAGE' ;
Thes = REDU dmthe_p (rvi . 'THE') ;
Thec = 'NOEL' $mthe_p (rvi . 'THE') ;
rhes muhes lhes cphes = PRPEAU prt ro_o mu_o cp_o lamb_o Thes ;
rhec muhec lhec cphec = PRPEAU prt ro_o mu_o cp_o lamb_o Thec ;
uhes = ELNO $mthe_p (KOPS qhe_p / rhec) ;
uhes = NOMC 'UZ' uhes ;
uhes = NOMC 'UZ' (-1 * uhes) ;
rvi . 'RCPHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' ('*' rhec cphec) ;
rvi . 'LHE_p' =  'KCHT' $mthe_p 'SCAL' 'CENTRE' lhec ;
rvi . 'UHE_p' =  'KCHT' $mthe_p 'VECT' 'SOMMET' uhes ;
 
 
fb1 = (exco uz (rvi . 'UHE')) ;
fbp1 = (exco uz (rvi . 'UHE_p')) ;
mess ;
mess ' -------------------------------' ;
mess 'rvi . uhe   : ' (mini fb1) (maxi fb1) ;
mess 'rvi . uhe_p : ' (mini fbp1) (maxi fbp1) ;
mess ' -------------------------------' ;
mess ;
 
 
MESS '  8) Calcul de la temperature dans l eau' ;
*   8) Calcul de la temperature dans l'eau
   EXEC rvhe ;
 
***Check de Convergence
 
   then = 'COPIER' (rvi . 'THE') ;
   tcbn = 'COPIER' (rvi . 'TCB') ;
   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 ;
****
 
'SI' ('EGA' maxiter vrai) ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'The computation has not converged') ;
   'ERREUR' 'Check calisa1.procedur !' ;
'FINSI' ;
 
********FIN DE LA BOUCLE DE CALCUL DIPHASIQUE de "essai11"**************
************************************************************************
 
tcpu9 ;
 
 
* CONTROLE DE LA TEMPERATURE DE PAROI : CAS DIPHASIQUE de "essai11"
************************************************************************
 
Tparoi  = 'NOEL' $mtif   (rvi . 'TCB') ;
chpobor 'Temperature de paroi' Tparoi ;
chpobor 'Temperature de saturation' Tpsat ;
errelT = ABS (KOPS (Tparoi - Tpsat) / Tpsat) ;
errelT = c_diph1 * errelT ;
 
*PRECISION A 1/1000 :
verifI = errelT MASQ 'INFERIEUR' 'SOMME' 0.001 ;
verifS = errelT MASQ 'SUPERIEUR' 'SOMME' 0.001 ;
mess verifS verifI ;
 
*PRECISION A 1/10000 :
verifI = errelT MASQ 'INFERIEUR' 'SOMME' 0.0001 ;
verifS = errelT MASQ 'SUPERIEUR' 'SOMME' 0.0001 ;
mess verifS verifI ;
 
*Erreur max :
errelTm = maxi errelT ;
mess errelTm ;
 
SI (errelTm > 1.e-4) ;
   MESS 'WARNING erreur sur la temperature paroi trop importante' ;
   MESS 'WARNING on stoppe le calcul' ;
 
SI interact ;
opti donn 5 ;
FINSI ;
 
FINSI ;
 
* FIN DU CONTROLE DE LA TEMPERATURE DE PAROI : CAS DIPHASIQUE "essai11"*
************************************************************************
 
tcpu9 ;
 
FINSI ; comm SI de la boucle monodiphasique ;
 
opti echo 1 ;
 
epsi0 = 1.e-6 ;
 
** Les Temperatures et les vitesses
 
tk = rvi.tcb ; comm temperatures solides en K ;
tc = tk - ck ; comm temperatures solides en C ;
t_tn = rvcb . 'INCO' . 'TCB' ;
 
Ts = (rvi . 'TCB') - CK ;
Tf = (rvi . 'THE') - CK ;
SpeedC = EXCO 'UZ' (rvi . 'UHE') ;
SpeedF = EXCO 'UZ' (rvi . 'UHE_p') ;
 
Tff = REDU dmthe_p Tf ;
Tfc = REDU dmthe Tf ;
 
opti echo 0 ;
 
* POST-TRAITEMENT et BILANS MONO et DIPHASIQUE
************************************************************************
 
* Calcul des densites de puissances aux interfaces fluide/solide
* --------------------------------------------------------------
 
** evfx sont des densites de puisances en W/m2 de type EVOLution
** F0x  sont des densites de puisances en W/m2 de type CHPOint
 
evf0 F0n = recal roug F0 $mtif mtif (0. (diamext / 2.)  z3 )
                     (0. (diamext / 2.)  z4 ) 50 ;
evf01 F01n = recal roug F01 $mtin_p mtin_p (0. (pani/ 2.)  z3 )
                     (0. (pani / 2.)  z4 ) 50 ;
evfn4 Fn4n = recal turq Fn4 $mtef_t mtef_t (rteexA 0. z3)
                    (rteexA 0. z4) 50 ;
 
IntF0  = maxi ( resu f0  ) ;
IntF01 = maxi ( resu f01 ) ;
IntF02 = maxi ( resu f02 ) ;
IntF03 = maxi ( resu f03 ) ;
IntF04 = maxi ( resu fn4 ) ;
 
phi0 = plin1 / ( pi * diamext ) ;
wtotcra = IntF0 ;
wpertpa = IntF01 ;
 
Lch = z4 - z3 ; comm longueur chauffante ;
 
* Pimp0 = puissance imposee sur la paroi seche du crayon
* ------------------------------------------------------
Pimpos0 = KOPS '*' ( DOMA $mtin 'XXVOLUM' ) pmsurf ;
Pimp0 = maxi (resu Pimpos0 ) ;
 
* Legendes dessins
* ----------------
tbb1 = TABLE;
tbb1.'TITRE' = TABLE ;
tbb1.1 = 'MARQ CROI REGU';
tbb1.'TITRE' . 1  = MOT 'Interne' ;
tbb1.2 = 'MARQ LOSA REGU';
tbb1.'TITRE' . 2  = MOT 'Externe' ;
 
 
tbb2 = TABLE;
tbb2.'TITRE' = TABLE ;
tbb2.1 = 'MARQ CROI REGU';
tbb2.'TITRE' . 1  = MOT 'TFI-Int' ;
tbb2.2 = 'MARQ LOSA REGU';
tbb2.'TITRE' . 2  = MOT 'TFI-Ext' ;
tbb2.3 = 'MARQ CARR REGU';
tbb2.'TITRE' . 3  = MOT 'TFE-Int' ;
tbb2.4 = 'MARQ ETOI REGU';
tbb2.'TITRE' . 4  = MOT 'TFE-Ext' ;
 
TCC1 = TABLE;
TCC1.'TITRE' = TABLE ;
TCC1.1 = 'MARQ CROI REGU';
TCC1.'TITRE' . 1  = MOT 'F. chaud' ;
TCC1.2 = 'MARQ LOSA REGU';
TCC1.'TITRE' . 2  = MOT 'F. froid' ;
 
 
* Temperatures tubes de force
* ---------------------------
 
opti elem seg2 ;
hori_E = DROI 15 ( rtiinA 0. 0. )  ( rteexA 0. 0. ) ;
evhor_E = evol chpo (INT_COMP (mt_tfi et mt_gaz et mt_tfe)
          tc hori_E ) hori_E ;
 
evgaii = (resex tc mtin)   coul roug ;
evgaie = (resex tc mtif)   coul rose ;
evtfii = (resex tc mtin_t) coul bleu  ;
evtfie = (resex tc mtif_t) coul vert ;
evtfei = (resex tc mten_t) coul roug ;
evtfee = (resex tc mtef_t) coul cyan ;
 
 
* Calcul des volumes
* ------------------
Volcb   =   DOMA $mtcb    'VOLUME' ;
Volcb_p =   DOMA $mtcb_p  'VOLUME' ;
Vol_tfi =   DOMA $mt_tfi  'VOLUME' ;
Vol_gaz =   DOMA $mt_gaz  'VOLUME' ;
Vol_tfe =   DOMA $mt_tfe  'VOLUME' ;
Volhe   =   DOMA $mthe    'VOLUME' ;
Volhe_p =   DOMA $mthe_p  'VOLUME' ;
VolcbAr =   DOMA $mtcbAr  'VOLUME' ;
VolcbEr =   DOMA $mtcbEr  'VOLUME' ;
Vol_lam =   DOMA $mt_lam  'VOLUME' ;
Vol_pe  =   DOMA $mt_pe   'VOLUME' ;
 
si egasol ;
 S_cbAr = SOMT ( KOPS VolcbAr * ( rvi . 'WgacbAr' ));
 S_cbEr = SOMT ( KOPS VolcbEr * ( rvi . 'WgacbEr' ));
 S_cb   = SOMT ( KOPS Volcb   * ( rvi . 'Wgacb'   ));
 S_cbp = SOMT ( KOPS Volcb_p * ( rvi . 'Wgacb_p' ));
 S_lam = SOMT ( KOPS Vol_lam * ( rvi . 'Wga_lam' ));
 S_pe  = SOMT ( KOPS Vol_pe  * ( rvi . 'Wga_pe'  ));
 S_tfi = SOMT ( KOPS Vol_tfi * ( rvi . 'Wga_tfi' ));
 S_gaz = SOMT ( KOPS Vol_gaz * ( rvi . 'Wga_gaz' ));
 S_tfe = SOMT ( KOPS Vol_tfe * ( rvi . 'Wga_tfe' ));
 S_tot1 =  S_cbAr + S_cbEr + S_cb + S_cbp + S_lam
      + S_pe + S_tfi + S_gaz + S_tfe ;
finsi ;
 
si egaliq ;
 S_he   = SOMT ( KOPS Volhe   * ( rvi . 'Wgahe'   ));
 S_hep = SOMT ( KOPS Volhe_p * ( rvi . 'Wgahe_p' ));
 S_tot2 =  S_he + S_hep ;
finsi ;
 
 
s_tot = 0. ;
si egasol ;
 s_tot = s_tot + S_tot1 ;
 sougams = 'Sources internes dans les structures' ;
sinon ;
 S_tot1 = 0. ;
 sougams = 'Sources internes dans les structures non compris ' ;
finsi ;
si egaliq ;
 S_tot = s_tot + S_tot2 ;
 sougaml = 'Sources internes dans les fluides' ;
sinon ;
 S_tot2 = 0. ;
 S_he = 0. ;
 S_hep = 0. ;
 sougaml = 'Sources internes dans les fluides non compris' ;
finsi ;
 
si ((non egasol) et (non egaliq)) ;
 sourgam = 'Sources internes non pris en compte' ;
sinon ;
 sourgam = 'Sources internes prises en compte' ;
finsi ;
 
 
 
* tff et tfc sont en degres CELSIUS
* ----------------------------------
*Fluide froid :
*--------------
hahep1 = mthe_p poin proc ( poin ( hahe_p chan poi1 ) 1 );
hahep2 = mthe_p poin proc ( poin (hahe_p  chan poi1 ) 2 );
hahep3 = mthe_p poin proc ( poin (hahe_p  chan poi1 ) 3 );
hahep4 = mthe_p poin proc ( poin (hahe_p  chan poi1 ) 4 );
tffe1 = ( maxi ( resu  ( redu tff hahep1 ) )) ;
tffe2 = ( maxi ( resu  ( redu tff hahep2 ) )) ;
tffe3 = ( maxi ( resu  ( redu tff hahep3 ) )) ;
tffe4 = ( maxi ( resu  ( redu tff hahep4 ) )) ;
tffe = ( tffe1 + tffe2 + tffe3 + tffe4  ) /
          (flottant ( nbno hahe_p )) ;
tffs0 = redu tff bahe_p ;
tffs = ( maxi ( resu tffs0 ) ) / (flottant ( nbno bahe_p )) ;
 
* Evolution fluide froid
* ----------------------
tff9 = redu tff mthe_p ;
tff9c = noel $mthe_p tff9 ;
mtff9c = doma $mthe_p centre ;
repe x (nbno mtff9c ) ;
pt = mtff9c point &x ;
val89 = coor 3 pt ;
wal89 = extr tff9c scal pt ;
 si ( ega &x 1 ) ;
  dqqf = val89 ;
  drrf = wal89 ;
 sinon ;
  dqqf = dqqf et val89 ;
  drrf = drrf et wal89 ;
 finsi ;
fin x ;
evtff9 = ( evol manu  'Z' dqqf drrf ) coul bleu ;
evtff9 = ordonner evtff9 ;
 
 
 
*Fluide chaud :
*--------------
tfce0 = redu tfc bahe ;
tfce = ( maxi ( resu tfce0 ) ) / (flottant ( nbno bahe )) ;
tfcs0 = redu tfc hahe ;
tfcs = ( maxi ( resu tfcs0 ) ) / (flottant ( nbno hahe )) ;
 
hffe = 'VARI' 'HLS' ppa ( tffe + ck ) ;
hffs = 'VARI' 'HLS' ppa ( tffs + ck ) ;
hfce = 'VARI' 'HLS' ppa ( tfce + ck ) ;
hfcs = 'VARI' 'HLS' ppa ( tfcs + ck ) ;
hlsat = 'VARI' 'HLS' ppa ;
 
* Evolution fluide chaud
* ----------------------
tfc9 = redu tfc mthe ;
tfc9c = noel $mthe tfc9 ;
mtfc9c = doma $mthe 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 ;
evtfc9 =  ( evol manu  'Z' dqqc drrc ) coul roug ;
evtfc9 = ordonner evtfc9 ;
 
* bilan global
*-------------
 
mDHhot  = (hfcs - hfce ) * qmhe  ;
mDHcold = (hffs - hffe ) * qmhe  ;
Win = Pimp0 + S_tot1 + S_tot2 - IntF04 ;
Wout = mDHcold + mDHhot ;
BilGlo = (S_tot + ( abs Pimp0) - IntF04) / Wout ;
 
 
* Calcul du flux critique a partir de la correlation de BOWRING
* -------------------------------------------------------------
 
tab1 = table ;
TAB1.'NIVEAU'   =0 ;
TAB1.'D_DIAM'   = Dhy ;  comm diam hydraulique fluide chaud ;
TAB1.'L_HEATED' = (Lch / 2.) ; comm demi longueur chauffante ;
TAB1.'P_IN'     = psatb * 1.e5 ;
TAB1.'G_IN'     = qhenom ; comm vitesse massique ;
TAB1.'T_IN'     = Theec ; comm temperature entree ;
bowr0 = bowring tab1 lat00 ;
bowr2 = bowr0 * pi * diamext * 0.01 ; comm W par cm ;
 
** Impression des resultats globaux
** --------------------------------
 
si ( non complet ) ;
 mess 'Calcul rapide monophasique seul' ;
 mess '*******************************' ;
 mess ;
finsi ;
 
ecform 'Pression                                         '
 '(F5.1)' 'bar'  psatb ;
ecform 'Diametre externe crayon chauffant                '
 '(F5.1)' 'mm'   ( 1000. * diamext) ;
ecform 'Longueur chauffante                              '
 '(F5.1)' 'mm'   ( 1000. * Lch ) ;
ecform 'Debit massique commun jambes chaude et froide    '
 '(F5.3)' 'kg/s' qmhe ;
ecform 'Puissance lineique maximale crayon chauffant     '
 '(F5.1)' 'W/cm' ( plin1 / 100. ) ;
ecform 'Flux critique (BOWRING) a mi-longueur chauffante '
 '(F5.1)' 'W/cm2' (bowr0/1.e4) ;
ecform 'Flux critique (BOWRING) a mi-longueur chauffante '
 '(F8.1)' 'W/cm' bowr2 ;
ecform 'Puissance incidente crayon  chauffant            '
 '(F8.1)' 'W'    Pimp0 ;
 
mess ;
mess sourgam ;
mess sougams ;
mess sougaml ;
si egasol ;
 ecform 'Sources internes structures                      '
  '(F8.1)' 'W' S_tot1  ;
finsi ;
si egaliq ;
 ecform 'Sources internes fluides                         '
  '(F8.1)' 'W' S_tot2 ;
finsi ;
ecform 'Sources internes totales                          '
  '(F8.1)' 'W' (S_tot1 + S_tot2) ;
mess ;
mess 'Puissance interface ';
ecform 'jambe chaude & parois separant les 2 ecoulements  '
  '(F8.1)' 'W' wpertpa ;
mess 'Puissance interface';
ecform 'jambe froide & parois separant les 2 ecoulements  '
  '(F8.1)' 'W' IntF02  ;
mess 'Puissance interface';
ecform 'jambe froide & tubes de force                     '
  '(F8.1)' 'W' IntF03  ;
ecform 'Pertes vers eau statique de la piscine            '
  '(F8.1)' 'W'  IntF04;
 
mess ;
Mess 'Fluide caloporteur : ';
mess '------------------' ;
ecform 'Temperature entree fluide chaud (imposee)            '
  '(F5.1)' 'degres C'  tfce ;
ecform 'Temperature sortie fluide chaud                      '
  '(F5.1)' 'degres C' tfcs ;
ecform 'Gradient de temperature fluide chaud                 '
  '(F5.1)' 'degres' (tfcs - tfce ) ;
 
ecform 'Temperature entree fluide froid                      '
  '(F5.1)' 'degres C' tffe  ;
 
ecform 'Temperature sortie fluide froid (imposee)            '
  '(F5.1)' 'degres C' tffs ;
 
ecform 'Gradient de temperature fluide froid                 '
  '(F5.1)' 'degres' (tffs - tffe )  ;
 
mess ;
Mess 'Temperature structures solides : ' ;
mess '------------------------------' ;
si complet ;
 ecform 'Temperature de saturation                            '
 '(F5.1)' 'degres C' (tsat - CK);
 ecform 'Temperature ebullition nucleee maximale (JENS & LOTTE) '
 '(F5.1)' 'degres C' (( maxi tpsat) - CK)  ;
finsi ; comm si complet ;
ecform 'Temperature maximale paroi externe crayon chauffant  '
 '(F5.1)' 'degres C' ( maxi ( redu tc mtif)) ;
 
ecform 'Temperature maximale paroi externe en contact piscine'
 '(F5.1)' 'degres C' ( maxi ( redu tc mtef_t)) ;
 
 
* TESTS RESULTATS et ERREURS
* --------------------------
vtest11 = -370.4 ; vtest12 = 581.1 ;
vtest21 = 1589.8 ; vtest22 = 638.4 ;
vtest31 = 2700.4 ; vtest32 = 2083.3 ;
vtest41 = 11600.8; vtest42 = 13150.1 ;
vtest51 =  252.0 ; vtest52 = 252.0 ;
vtest61 =  269.5 ; vtest62 = 270.6 ;
vtest71 =  245.5 ; vtest72 = 247.2  ;
vtest81 =  251.0 ; vtest82 = 251.0 ;
 
si complet ;
 ecrela1 = ABS ((Wpertpa - vtest11) / vtest11) ;
 ecrela2 = ABS ((IntF02  - vtest21) / vtest21) ;
 ecrela3 = ABS ((IntF03  - vtest31) / vtest31) ;
 ecrela4 = ABS ((IntF04  - vtest41) / vtest41) ;
 ecrela5 = ABS ((tfce    - vtest51) / vtest51) ;
 ecrela6 = ABS ((tfcs    - vtest61) / vtest61) ;
 ecrela7 = ABS ((tffe    - vtest71) / vtest71) ;
 ecrela8 = ABS ((tffs    - vtest81) / vtest81) ;
sinon ;
 ecrela1 = ABS ((Wpertpa - vtest12) / vtest12) ;
 ecrela2 = ABS ((IntF02  - vtest22) / vtest22) ;
 ecrela3 = ABS ((IntF03  - vtest32) / vtest32) ;
 ecrela4 = ABS ((IntF04  - vtest42) / vtest42) ;
 ecrela5 = ABS ((tfce    - vtest52) / vtest52) ;
 ecrela6 = ABS ((tfcs    - vtest62) / vtest62) ;
 ecrela7 = ABS ((tffe    - vtest72) / vtest72) ;
 ecrela8 = ABS ((tffs    - vtest82) / vtest82) ;
finsi ;
 
epseca = 1D-1 ;
testsym1 = > ecrela1 epseca ;
testsym2 = > ecrela2 epseca ;
testsym3 = > ecrela3 epseca ;
testsym4 = > ecrela4 epseca ;
testsym5 = > ecrela5 epseca ;
testsym6 = > ecrela6 epseca ;
testsym7 = > ecrela7 epseca ;
testsym8 = > ecrela8 epseca ;
 
testerr = testsym1 'OU' testsym2 'OU' testsym3 'OU' testsym4
   'OU' testsym5 'OU' testsym6 'OU' testsym7 'OU' testsym8 ;
 
'SI' testerr ;
   'MESSAGE' 'Il y a eu des erreurs' ;
   'ERREUR' 5 ;
'FINSI' ;
 
 
SI graph ;
 
opti echo 1 ;
 
si graffic ;
  'OPTION' 'TRAC' 'PSC' ;
*  'OPTION' 'ECHO' 0 ;
'SINON' ;
  'OPTION' 'TRAC' 'X' ;
'FINSI' ;
 
 
 loupe mtif (t_tn - CK) 10.
      'Temperature externe crayon chauffant (C)' SURF 0  ;
 
 dess evf0 mima titx 'cote axiale m'
 tity 'Flux (W/m2)'
 titr 'Interface fluide chaud & crayon chauffant' ;
 
 dess evf01 mima titx 'cote axiale m'
 tity 'Flux (W/m2)'
 titr 'Interface fluide chaud & parois separant les 2 ecoulements' ;
 
  dess evfn4 mima titx 'cote axiale m'
  tity 'Flux (W/m2)'
  titr 'Interface tube de force & piscine' ;
 dess ( evtfc9 et evtff9 ) mima titx 'm' tity 'degres C'
 titre 'Temperatures fluides' lege TCC1 ;
 
 dess (evgaii et evgaie) mima titx 'm'
 tity 'Degres C' titr 'Temperatures gaine' LEGE tbb1 ;
 
 dessin (evtfii et evtfie et evtfei et evtfee) mima
       titx 'm' tity 'Degres C'
       titre 'Temperatures tube de force interne sur axe vertical'
       LEGE tbb2 ;
 dess evhor_E titx 'm' tity 'Degres C'
 titr 'Temperatures tubes de force sur axe horizontal a z=0';
 
 
FINSI ; comm si graph ;
 
****************FIN DU POST-TRAITEMENT**********************************
************************************************************************
 
 
 
si ( non interact ) ;
 
chpvide matvide = 'KOPS' 'MATRIK' ;
rvcb . METHINV.'MATASS'=matvide;
rvcb . METHINV.'MAPREC'=matvide;
rvhe . METHINV.'MATASS'=matvide;
rvhe . METHINV.'MAPREC'=matvide;
menage;
* En principe il n'y a plus de matrik !!!
 fin ;
finsi ;