Télécharger la plaquette Cast3M


fr en

Télécharger pq1xx.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  pq1.dgibi
  2. *---------------------------------------------------------
  3. * Cas test de non régression sur le controle des maillages
  4. *---------------------------------------------------------
  5. *     Pressurisation d'une enceinte type Phébus
  6. *
  7. * Le maillage correspond à une enceinte cylindrique
  8. * d'environ 10 m3 avec un mur en contact avec la
  9. * paroi verticale (10 cm d'acier)
  10. * Tout le volume est initialement a 1bar et 40oC
  11. * et la température du mur est mise à 60oC
  12. * au début du calcul. On calcule la pressurisation
  13. * de cette enceinte sur 50 secondes en injectant un débit
  14. * de 50g/s de vapeur à 150oC. Ce test (un peu long)
  15. * verifie le demarrage de la condensation a 20"
  16. * verifie l'evolution moyenne de la temprature gaz
  17. * verifie la pression max a 50"
  18. * verifie la vitesse max a 50" (Convection naturelle)
  19. * verifie la masse d'eau a 50"
  20. * Auteurs : E. Studer,                 Novembre 1999
  21. *                             revisite Mars     2001
  22. *--------------------------------------------------------
  23. COMPLET = FAUX ;
  24. *COMPLET = VRAI ;
  25.  GRAPH =  VRAI                                          ;
  26.  GRAPH =  FAUX                                          ;
  27.  'OPTI' 'TRAC' 'X'                                ;
  28. *'OPTI' 'TRAC' 'PSC'                              ;
  29.  
  30.  'SI' COMPLET ;
  31.  
  32.  nbit=100;
  33.  DT0 = 1. ;
  34.  n1   =  1          ; n2   =  4          ; n3   =  4          ;
  35.  n4   =  8          ; nn   =  2          ;
  36.  
  37.  'SINON' ;
  38.  
  39.  nbit= 5 ;
  40.  DT0 = 10. ;
  41.  n1   =  1          ; n2   =  2          ; n3   =  4          ;
  42.  n4   =  4          ; nn   =  1          ;
  43.  
  44.  'FINSI' ;
  45.  
  46. *--------------------------------------------------------
  47. * Definition du maillage de l'enceinte cylindrique
  48. *
  49.  'OPTI' 'DIME' 3 'ELEM' 'CU20'                                ;
  50.  
  51.  ri   = 1.052      ; sp   =  0.10     ; hc   = 4.163      ;
  52.  
  53.  
  54.  
  55.  epsi =  1.000e-2   ;                                         ;
  56.  
  57.  p0   =     0.000             0.000                0.000      ;
  58.  px   = -1000.000             0.000                0.000      ;
  59.  py   =     0.000         -1000.000                0.000      ;
  60.  pz   =     0.000             0.000             1000.000      ;
  61.  cd   =     0.000             0.000              -20.000      ;
  62.  ph0  =     0.000             0.000                hc         ;
  63.  phx  =     ri                0.000                hc         ;
  64.  phy  =     0.000             ri                   hc         ;
  65.  
  66.  fg1 = 0.25  ;
  67.  fg2 = fg1 * (2.0 ** 0.5) / 2. ;
  68.  
  69.  p1   =  (ri*fg1)             0.000                0.000      ;
  70.  p2   =  (ri*fg2)             (ri*fg2)             0.000      ;
  71.  p3   =  0.000                (ri*fg1)             0.000      ;
  72.  p4   =   ri                  0.000                0.000      ;
  73.  p5   =  0.000                 ri                  0.000      ;
  74.  p6   =  (ri+sp)              0.000                0.000      ;
  75.  p7   =  0.000                (ri+sp)              0.000      ;
  76.  
  77. * Hauteur de l'enceinte
  78.   h1 = 4.163                                                    ;
  79. * Vecteur de translation
  80.   v1 = 0. 0. h1                                               ;
  81.  
  82.  l1   =  'DROI' p0 p1 n1                                        ;
  83.  l2   =  'DROI' p1 p2 n1                                        ;
  84.  l3   =  'DROI' p2 p3 n1                                        ;
  85.  l4   =  'DROI' p3 p0 n1                                        ;
  86.  l5   =  'CERC' p4 p0 p5 (2*n1)                                 ;
  87.  l6   =  'CERC' p6 p0 p7 (2*n1)                                 ;
  88.  
  89.  basf0=  'DALL' l1 l2 l3 l4 'PLAN'                              ;
  90.  basf1=('REGL' (l2 'ET' l3) l5 n2) ;
  91.  l44= cote 2 basf1;
  92.  ax4= (inve l4) et l44 ;
  93.  l11= cote 4 basf1;
  94.  ax1= l11 et (inve l1) ;
  95.  
  96.  basf =  basf0 'ET' ('REGL' (l2 'ET' l3) l5 n2) ;
  97.  'ELIM' basf epsi ;
  98.  basf =  basf 'ET' ('SYME'  basf 'DROI' p0 p3) ;
  99.  ax11 =  ('SYME'  ax1  'DROI' p0 p3) 'ET' (inve ax1) ;
  100.  'ELIM' basf epsi ;
  101.  basf =  basf 'ET' ('SYME'  basf 'DROI' p0 p1) ;
  102.  ax44 =  (inve ax4)  'ET' ('SYME'  ax4  'DROI' p0 p4) ;
  103.  'ELIM' basf epsi ;
  104.  
  105.  basm =  'REGL' l5 l6 n3                                        ;
  106.  basm =  basm 'ET' ('SYME'  basm 'DROI' p0 p3) ;
  107.  'ELIM' basm epsi ;
  108.  basm =  basm 'ET' ('SYME'  basm 'DROI' p0 p1) ;
  109.  'ELIM' basm epsi ;
  110.  
  111. * Creation du volume
  112.  
  113.  dx = ri / 2. ;
  114.  nz1 = ('ENTIER' ( h1 / dx ))*nn ;
  115.  
  116.  bas = basf 'VOLU' nz1 'TRAN' v1 ;
  117.  mbas = basm 'VOLU' nz1 'TRAN' v1 ;
  118.  mbas = mbas 'COUL' 'ROUG' ;
  119.  plan1 = ax11 'TRAN' nz1 v1 ;
  120.  plan4 = ax44 'TRAN' nz1 v1 ;
  121.  'ELIM' (bas et plan1 et plan4) epsi ;
  122.  
  123.  mt = bas ;
  124.  wall = mbas ;
  125.  elim (mt et wall) epsi ;
  126.  
  127. * Localisation d'une brèche éventuelle au bas de l'enceinte
  128.  
  129.  pjg  =  'POIN' basf 'PROC' (0.000    0.000    0.000)             ;
  130.  jg   = ('ELEM' basf 'APPUIE' 'LARGEMENT' pjg) 'COUL' 'VERT'      ;
  131.  
  132. *--------------------------------------------------------------------
  133. * Fin de la définition du maillage
  134. *--------------------------------------------------------------------
  135.  
  136.  
  137. *--------------------------------------------------------------------
  138. * Début de l'initialisation de la procédure ENCEINTE : table RXT
  139. *--------------------------------------------------------------------
  140.  
  141.  rxt             = 'TABLE' ;
  142.  rxt . 'VERSION' = 'V0' ;
  143.  
  144. *-- Nom du volume fluide
  145.  rxt . 'vtf'     = mt ;
  146.  rxt . 'epsi'    = epsi ;
  147.  
  148. *-- Definition des murs de l'enceinte : ici un seul mur
  149. *-- en ACIER dont on traitera la thermique dans l'épaisseur
  150. *-- et que l'on initialise a 40oC
  151.  
  152. *-- On definit d'abord le matériau ACIER avec sa conductivite
  153. *-- thermique LAMBDA (W/m/K) et le produit ro*Cp (J/m3/K)
  154.  rxt . 'THERMP' = VRAI                                        ;
  155.  rxt . 'vtp'    = wall                                        ;
  156.  rxt . 'LAMBDA' = 15.                                         ;
  157.  rxt . 'ROCP'   = 3.9E6                                       ;
  158.  rxt . 'Tp0'    = 60.                                         ;
  159.  rxt . 'ECHAN'  = 10.                                         ;
  160.  
  161. *-- On definit un point interne au maillage pour imposer la valeur de
  162. *-- la pression
  163.  rxt . 'pi'   =  (0.0 0.0 0.5)         ;
  164.  
  165. *-- On indique que le calcul comporte de la vapeur d'eau
  166.  rxt . 'VAPEUR' = VRAI                                        ;
  167.  rxt . 'Yvap0'  = 0.0023                      ;
  168.  
  169. *-- Conditions initiales dans l'enceinte de test
  170.  rxt . 'TF0'    = 40.0                        ;
  171.  rxt . 'PT0'    = 1.00000e5                   ;
  172.  
  173. *-- Definition du scenario à la breche
  174.  
  175.  rxt . 'Breches' = 'TABLE' ;
  176.  rxt . 'Breches' . 'A' = 'TABLE' ;
  177.  rxt . 'Breches' . 'A' . 'Maillage' =  jg ;
  178.  rxt . 'Breches' . 'A' . 'diru'     =  (0. 0. 1.) ;
  179.  rxt . 'Breches' . 'A' . 'scenario' = 'TABLE' ;
  180.  rxt . 'Breches' . 'A' . 'scenario' . 't'    =  PROG 0.0  1000. ;
  181.  rxt . 'Breches' . 'A' . 'scenario' . 'qair' =  PROG 0.0  0.0   ;
  182.  rxt . 'Breches' . 'A' . 'scenario' . 'qeau' =  PROG 0.05 0.05  ;
  183.  rxt . 'Breches' . 'A' . 'scenario' . 'tinj' =  PROG 150. 150.  ;
  184.  
  185. *-- On impose la viscosite turbulente (m2/s)
  186.  rxt . 'MODTURB' = 'NUTURB' ;
  187.  rxt . 'NUT'     =  1.e-2   ;
  188.  
  189. *-- On impose le pas de temps (s)
  190.  rxt . 'DT0'    =  DT0                                          ;
  191.  
  192. *-- On active le recalcul automatique du préconditionnement
  193. *-- toutes les 5 itérations
  194.  rxt . 'FRPREC' = 5    ;
  195.  rxt . 'DETMAT' = VRAI ;
  196. *rxt.'TRTF'   = FAUX                                       ;
  197.  
  198.  
  199. Si VRAI;
  200. *-- On lance le calcul sur 20 itérations d'une seconde
  201.  rxt . 'GRAPH' = GRAPH ;
  202.  
  203.  EXECRXT  2    rxt                              ;
  204.  EXECRXT  (nbit - 2) rxt                        ;
  205.  
  206. list rxt.'TIC'.'Tfm' ;
  207. list rxt.'TIC'.'PT'  ;
  208. list rxt.'TIC'.'Qc'  ;
  209. list rxt.'TIC'.'LMAXU';
  210.  
  211. ltfm = 'PROG' 
  212. 40.000    65.434    73.712    81.158    86.809    90.721  ;
  213. lPT  = 'PROG' 
  214. 1.00000E+05 1.06009E+05 1.21199E+05 1.29592E+05 1.36461E+05 1.43987E+05;
  215. Lqc  = 'PROG' 
  216. 0.0000    0.0000    0.0000   4.06133E-04 2.65247E-03 4.88852E-03;
  217. Lmaxu = 'PROG' 
  218. 0.0000   0.81320    2.0813    2.6901    2.2350    2.5277  ;
  219.  
  220. tic=rxt.'TIC' ;
  221. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  222. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  223. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  224. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  225.  
  226. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  227.  
  228. ierr = 0 ;
  229. Si (ERtf '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  230. Si (ERPT '>' 1.e-3) ; ierr = ierr + 1 ; Finsi ;
  231. Si (ERQc '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  232. Si (ERum '>' 1.e-2) ; ierr = ierr + 1 ; Finsi ;
  233.  
  234. ******************************************************************
  235.  
  236.  rxt.'vtf'  =  (chan 'MACRO' mt)           ;
  237.  rxt.'vtp'  =  (chan 'MACRO' wall)          ;
  238.  rxt . 'Breches' . 'A' . 'Maillage' = (chan 'MACRO' jg) ;
  239.  
  240.  rxt = enlev rxt 'GEO';
  241.  rxt = enlev rxt 'TBT';
  242.  rxt = enlev rxt 'TIC';
  243.  
  244.  EXECRXT 2  rxt               ;
  245.  EXECRXT (nbit - 2) rxt            ;
  246.  
  247. list rxt.'TIC'.'Tfm' ;
  248. list rxt.'TIC'.'PT' ;
  249. list rxt.'TIC'.'Qc' ;
  250. list rxt.'TIC'.'LMAXU';
  251.  
  252. tic=rxt.'TIC' ;
  253. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  254. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  255. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  256. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  257.  
  258. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  259.  
  260. Si (ERtf '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  261. Si (ERPT '>' 1.e-3) ; ierr = ierr + 1 ; Finsi ;
  262. Si (ERQc '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  263. Si (ERum '>' 1.e-2) ; ierr = ierr + 1 ; Finsi ;
  264.  
  265.  
  266. ******************************************************************
  267.  
  268.  rxt.'vtf'  =  (chan 'QUAF' mt)             ;
  269.  rxt.'vtp'  =  (chan 'QUAF' wall)            ;
  270.  rxt . 'Breches' . 'A' . 'Maillage' = (chan 'QUAF' jg) ;
  271.  rxt = enlev rxt 'GEO';
  272.  rxt = enlev rxt 'TBT';
  273.  rxt = enlev rxt 'TIC';
  274.  
  275.  EXECRXT 2  rxt               ;
  276.  EXECRXT (nbit - 2) rxt            ;
  277.  
  278. list rxt.'TIC'.'Tfm' ;
  279. list rxt.'TIC'.'PT' ;
  280. list rxt.'TIC'.'Qc' ;
  281. list rxt.'TIC'.'LMAXU';
  282.  
  283.  
  284. tic=rxt.'TIC' ;
  285. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  286. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  287. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  288. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  289.  
  290. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  291.  
  292. Si (ERtf '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  293. Si (ERPT '>' 1.e-3) ; ierr = ierr + 1 ; Finsi ;
  294. Si (ERQc '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  295. Si (ERum '>' 1.e-2) ; ierr = ierr + 1 ; Finsi ;
  296.  
  297. finsi;
  298.  
  299. ******************************************************************
  300.  rxt.'vtf'  =  (chan 'LINE' mt)             ;
  301.  rxt.'vtp'  =  (chan 'LINE' wall)            ;
  302.  rxt . 'Breches' . 'A' . 'Maillage' = (chan 'LINE' jg) ;
  303.  rxt = enlev rxt 'GEO';
  304.  rxt = enlev rxt 'TBT';
  305.  rxt = enlev rxt 'TIC';
  306.  
  307.  EXECRXT 2  rxt               ;
  308.  EXECRXT (nbit - 2) rxt            ;
  309.  
  310. list rxt.'TIC'.'Tfm' ;
  311. list rxt.'TIC'.'PT' ;
  312. list rxt.'TIC'.'Qc' ;
  313. list rxt.'TIC'.'LMAXU';
  314.  
  315. ltfm = 'PROG' 
  316. 40.000    67.363    75.733    83.257    88.699    92.375  ;
  317. lPT  = 'PROG' 
  318. 1.00000E+05 1.06452E+05 1.22785E+05 1.31710E+05 1.39047E+05 1.46818E+05;
  319. Lqc  = 'PROG' 
  320. 0.0000    0.0000    0.0000   1.04818E-03 3.74210E-03 6.39247E-03;
  321. Lmaxu = 'PROG' 
  322. 0.0000   0.81320    2.0604    2.5919    2.1446    2.4526  ;
  323.  
  324. tic=rxt.'TIC' ;
  325. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  326. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  327. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  328. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  329.  
  330. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  331.  
  332. Si (ERtf '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  333. Si (ERPT '>' 1.e-3) ; ierr = ierr + 1 ; Finsi ;
  334. Si (ERQc '>' 1.e-4) ; ierr = ierr + 1 ; Finsi ;
  335. Si (ERum '>' 1.e-2) ; ierr = ierr + 1 ; Finsi ;
  336.  
  337.  
  338.  
  339. Si GRAPH ;
  340.  
  341. $vtf=rxt.'GEO'.'$vtf' ;
  342. mt =doma $vtf maillage ;
  343. Mpl1=chan 'QUAF' plan1 ;
  344. Mpl4=chan 'QUAF' plan4 ;
  345. ELIM (mt et Mpl1 et Mpl4) epsi ;
  346. $mpl1= mode Mpl1 'NAVIER_STOKES' MACRO ;
  347. $mpl4= mode Mpl4 'NAVIER_STOKES' MACRO ;
  348. plan1= doma $mpl1 maillage ;
  349. plan4= doma $mpl4 maillage ;
  350. plan=plan1 et plan4 ;
  351.  
  352. paroif = rxt.'GEO'.'paroif';
  353. rho=rxt.'TIC'.'RHO';
  354. rvp=rxt.'TIC'.'RVP';
  355. tf=rxt.'TIC'.'TF';
  356. un=rxt.'TIC'.'UN' ;
  357. un1=redu un plan ;
  358. ung= vect un1 1. ux uy uz jaune ;
  359. trace ung plan ;
  360. opti isov suli ;
  361. trace rho plan 'TITRE'' Rho' ;
  362. trace rvp plan 'TITRE'' Rvp' ;
  363. trace tf plan 'TITRE'' Tf ' ;
  364.  
  365. trace rho paroif 'TITRE'' Rho' ;
  366. trace rvp paroif 'TITRE'' Rvp' ;
  367. trace Tf paroif 'TITRE'' Tf ' ;
  368. Fcond = tic.'Fcondw';
  369. trace Fcond paroif 'TITRE' ' Flux de condensation Kg / m**2 ' ;
  370.  
  371.  axe = p0 d nz1 (p0 plus v1) ;
  372.  axe = chan axe 'QUAF' ;
  373.  elim (axe et mt) epsi ;
  374.  
  375. evr= evol 'CHPO' rho axe ;
  376. dess evr 'TITRE' ' Rho axe ';
  377.  
  378. evt = evol 'CHPO' Tf axe ;
  379. dess evt 'TITRE' ' Température axe ';
  380.  
  381. evvp= evol 'CHPO' rvp axe ;
  382. dess evvp 'TITRE' ' rvp axe ';
  383.  
  384.  
  385.  
  386. 'FINSI' ;
  387.  
  388.  
  389. 'SI' (IERR '>' 0) ;
  390.  'MESS' 'Il y a des problemes !!!' ;
  391.  'ERRE' 5 ;
  392. 'SINO';
  393.  'MESS' 'Tout s est bien passe!' ;
  394. 'FINS' ;
  395.  
  396.  
  397.  'FIN' ;
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales