Accueil | Cast3M

* fichier :  pressutq.dgibi
*---------------------------------------------------------
*     Pressurisation d'une enceinte type Phébus
*
* Le maillage correspond à une enceinte cylindrique
* d'environ 10 m3 avec un mur en contact avec la
* paroi verticale (10 cm d'acier)
* Tout le volume est initialement a 1.3bar et 123oC
* et la température du mur est mise à 100oC
* au début du calcul. On calcule la pressurisation
* de cette enceinte sur 50 secondes en injectant un débit
* de 10g/s de vapeur à 130oC. Ce test (un peu long)
* verifie le demarrage de la condensation a 20"
* verifie l'evolution moyenne de la temprature gaz
* verifie la pression max a 50"
* verifie la vitesse max a 50" (Convection naturelle)
* Auteurs : E. Studer, J.P. Magnaud    Novembre 1999
*                             revisite Avril    2002
*--------------------------------------------------------
COMPLET = FAUX ;
*COMPLET = VRAI ;
 GRAPH =  VRAI                                          ;
 GRAPH =  FAUX                                          ;
 
 'SI' COMPLET ;
 
 nbit=100;
 DT0 = 1. ;
 n1   =  1          ; n2   =  4          ; n3   =  4          ;
 n4   =  8          ; nn   =  2          ;
 
 'SINON' ;
 
 nbit= 5 ;
 DT0 = 10. ;
 n1   =  1          ; n2   =  2          ; n3   =  4          ;
 n4   =  4          ; nn   =  1          ;
 
 'FINSI' ;
 
*--------------------------------------------------------
* Definition du maillage de l'enceinte cylindrique
*
 'OPTI' 'DIME' 3 'ELEM' 'CU20'                      ;
 
 ri   = 1.052      ; sp   =  0.10     ; hc   = 4.163      ;
 
 
 
 epsi =  1.000e-2   ;                                         ;
 
 p0   =     0.000             0.000                0.000      ;
 px   = -1000.000             0.000                0.000      ;
 py   =     0.000         -1000.000                0.000      ;
 pz   =     0.000             0.000             1000.000      ;
 cd   =     0.000             0.000              -20.000      ;
 ph0  =     0.000             0.000                hc         ;
 phx  =     ri                0.000                hc         ;
 phy  =     0.000             ri                   hc         ;
 
 fg1 = 0.25  ;
 fg2 = fg1 * (2.0 ** 0.5) / 2. ;
 
 p1   =  (ri*fg1)             0.000                0.000      ;
 p2   =  (ri*fg2)             (ri*fg2)             0.000      ;
 p3   =  0.000                (ri*fg1)             0.000      ;
 p4   =   ri                  0.000                0.000      ;
 p5   =  0.000                 ri                  0.000      ;
 p6   =  (ri+sp)              0.000                0.000      ;
 p7   =  0.000                (ri+sp)              0.000      ;
 
* Hauteur de l'enceinte
  h1 = 4.163                                                    ;
* Vecteur de translation
  v1 = 0. 0. h1                                               ;
 
 l1   =  'DROI' p0 p1 n1                                        ;
 l2   =  'DROI' p1 p2 n1                                        ;
 l3   =  'DROI' p2 p3 n1                                        ;
 l4   =  'DROI' p3 p0 n1                                        ;
 l5   =  'CERC' p4 p0 p5 (2*n1)                                 ;
 l6   =  'CERC' p6 p0 p7 (2*n1)                                 ;
 
 basf0=  'DALL' l1 l2 l3 l4 'PLAN'                              ;
 basf1=('REGL' (l2 'ET' l3) l5 n2) ;
 l44= cote 2 basf1;
 ax4= (inve l4) et l44 ;
 l11= cote 4 basf1;
 ax1= l11 et (inve l1) ;
 
 basf =  basf0 'ET' ('REGL' (l2 'ET' l3) l5 n2) ;
 'ELIM' basf epsi ;
 basf =  basf 'ET' ('SYME'  basf 'DROI' p0 p3) ;
 ax11 =  ('SYME'  ax1  'DROI' p0 p3) 'ET' (inve ax1) ;
 'ELIM' basf epsi ;
 basf =  basf 'ET' ('SYME'  basf 'DROI' p0 p1) ;
 ax44 =  (inve ax4)  'ET' ('SYME'  ax4  'DROI' p0 p4) ;
 'ELIM' basf epsi ;
 
 basm =  'REGL' l5 l6 n3                                        ;
 basm =  basm 'ET' ('SYME'  basm 'DROI' p0 p3) ;
 'ELIM' basm epsi ;
 basm =  basm 'ET' ('SYME'  basm 'DROI' p0 p1) ;
 'ELIM' basm epsi ;
 
* Creation du volume
 
 dx = ri / 2. ;
 nz1 = ('ENTIER' ( h1 / dx ))*nn ;
 
 bas = basf 'VOLU' nz1 'TRAN' v1 ;
 mbas = basm 'VOLU' nz1 'TRAN' v1 ;
 mbas = mbas 'COUL' 'ROUG' ;
 plan1 = ax11 'TRAN' nz1 v1 ;
 plan4 = ax44 'TRAN' nz1 v1 ;
 'ELIM' (bas et plan1 et plan4) epsi ;
 
 mt = bas ;
 wall = enve bas ;
 mur  = DIFF wall basf;
 mur =mur 'COUL' 'ROUG' ;
 
 elim (mt et mur) epsi ;
 
 
* Localisation d'une brèche éventuelle au bas de l'enceinte
 
 pjg  =  'POIN' basf 'PROC' (0.000    0.000    0.000)             ;
 jg   = ('ELEM' basf 'APPUIE' 'LARGEMENT' pjg) 'COUL' 'VERT'      ;
 
*--------------------------------------------------------------------
* Fin de la définition du maillage
*--------------------------------------------------------------------
 
*--------------------------------------------------------------------
* Début de l'initialisation de la procédure ENCEINTE : table RXT
*--------------------------------------------------------------------
 
 rxt          =  'TABLE'                                      ;
 
*-- Nom du volume fluide
 rxt.'vtf'    =   mt                                        ;
 rxt.'epsi'   =   epsi                                      ;
 
* Definition du bidon :
 rxt.'TIMP1'             = table                             ;
 rxt.'TIMP1'.'MAILLAGE'  = mur                               ;
 rxt.'TIMP1'.'t'         = prog   0.0    3000.               ;
 rxt.'TIMP1'.'TIMP'      = prog   100.0  100.0               ;
 rxt.'TIMP1'.'ECHAN'     = 10.                               ;
 
*-- Conditions initiales dans l'enceinte de test
 rxt.'TF0'    =      123.0                         ;
 rxt.'PT0'    =        1.30000e5                   ;
 rxt.'Yvap0'  =        1.0E-8                      ;
 
*-- On positionne une brèche
 rxt.'breche' =   jg                                         ;
 rxt.'diru1'  =   0. 0. 1.                                   ;
 
*-- On definit un point interne au maillage pour imposer la valeur de
*-- la pression
 rxt.'pi'     =  (0.0 0.0 0.5)                               ;
 
*-- On indique que le calcul comporte de la vapeur d'eau
 rxt.'VAPEUR' =  VRAI                                        ;
 
*-- On active le recalcul automatique du préconditionnement
*-- toutes les 5 itérations
 rxt.'FRPREC' = 5                                          ;
 rxt.'DETMAT' = VRAI                                       ;
 
*-- Definition du scenario thermohydraulique
 rxt.'scenario'          = table                              ;
 
*-- Conditions a la breche (Obligatoire pour l'instant)
 rxt.'scenario'.'t'      = prog   0.0 1000.0                  ;
 rxt.'scenario'.'qeau'   = prog   0.010  0.010                ;
 rxt.'scenario'.'qair'   = prog   0.000  0.000                ;
 rxt.'scenario'.'tinj'   = prog   130.0  130.0                ;
 
*-- On impose le pas de temps (s)
 rxt.'DT0'    =  DT0                                          ;
 
*-- On impose la viscosite turbulente (m2/s)
 rxt.'MODTURB'='NUTURB'                             ;
 rxt.'NUT'    =  1.e-2                              ;
 
*-- On lance le calcul sur 20 itérations d'une seconde
 rxt.'GRAPH'=GRAPH                      ;
 
 execrxt 2    rxt                              ;
 execrxt (nbit - 2) rxt                        ;
 
 
list rxt.TIC.'Tfm' ;
list rxt.TIC.'PT'  ;
list rxt.TIC.'Qc'  ;
list rxt.TIC.'LMAXU';
 
 'SI' (NON  COMPLET) ;
 
ltfm=Prog 123.00       120.38       122.24       118.51
          116.67       115.06;
lPT =Prog 1.30000E+05  1.30840E+05  1.30301E+05  1.32429E+05
          1.32769E+05  1.33224E+05;
Lqc =Prog 0.0000       0.0000       0.0000       0.0000
          0.0000       0.0000;
Lmaxu=Prog 0.0000      0.11919      0.12541       1.5872
           1.0398       1.9797;
 
tic=rxt.'TIC' ;
ERtf=somm( abs (ltfm - tic.'Tfm') )/ 110. ;
ERPT=somm( abs (lPT  - tic.'PT' ) ) /1.e5 ;
ERQc=somm( abs (lqc  - tic.'Qc' ) )  ;
ERum=somm( abs (Lmaxu  - tic.'LMAXU' ) )/2. ;
 
Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
 
Si (ERtf '>' 1.e-3) ; erreur 5 ; Finsi ;
Si (ERPT '>' 1.e-4) ; erreur 5 ; Finsi ;
Si (ERQc '>' 1.e-4) ; erreur 5 ; Finsi ;
Si (ERum '>' 1.e-2) ; erreur 5 ; Finsi ;
 
 'FINSI'  ;
 
Si GRAPH ;
 
$vtf=rxt.'GEO'.'$vtf' ;
mt  =doma $vtf maillage ;
Mpl1=chan 'QUAF' plan1 ;
Mpl4=chan 'QUAF' plan4 ;
ELIM (mt et Mpl1 et Mpl4) epsi ;
$mpl1= mode Mpl1 'NAVIER_STOKES' MACRO ;
$mpl4= mode Mpl4 'NAVIER_STOKES' MACRO ;
plan1= doma $mpl1 maillage ;
plan4= doma $mpl4 maillage ;
plan=plan1 et plan4 ;
 
paroif = rxt.'GEO'.'mtp1';
rho=rxt.'TIC'.'RHO';
rvp=rxt.'TIC'.'RVP';
tf=rxt.'TIC'.'TF';
un=rxt.'TIC'.'UN' ;
un1=redu un plan ;
ung= vect un1 1. ux uy uz jaune ;
trace ung plan ;
opti isov suli ;
trace rho plan 'TITRE'' Rho' ;
trace rvp plan 'TITRE'' Rvp' ;
trace tf  plan 'TITRE'' Tf ' ;
 
trace rho paroif  'TITRE'' Rho' ;
trace rvp paroif  'TITRE'' Rvp' ;
trace Tf  paroif  'TITRE'' Tf ' ;
Fcond = rxt.'TIC'.'Fcond1';
trace Fcond paroif  'TITRE' ' Flux de condensation Kg / m**2 ' ;
 
 axe  = p0 d nz1 (p0 plus v1) ;
 axe = chan axe 'QUAF' ;
 elim (axe et mt) epsi ;
 
evr= evol 'CHPO' rho axe ;
dess evr 'TITRE' ' Rho axe ';
 
evt = evol 'CHPO' Tf  axe ;
dess evt 'TITRE' ' Température axe ';
 
evvp= evol 'CHPO' rvp axe ;
dess evvp 'TITRE' ' rvp axe ';
 
'FINSI' ;
 
 'FIN' ;