Télécharger la plaquette Cast3M


fr en

Télécharger carre_therper.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  carre_therper.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ******************************************************************
  5. * CALCUL DU QUARRE A CHOC                                        *
  6. * Modele: gaz multi-especes "thermally perfect"                  *
  7. * Gaz multi-especes "thermally   perfect"                        *
  8. * Controle de la symetrie                                        *
  9. *                                                                *
  10. * A. BECCANTINI,    SEMT/LTMF,     MARS 1999                     *
  11. ******************************************************************
  12.  
  13.  
  14. 'OPTION'  'DIME' 2 'ELEM' 'QUA4' 'ISOV' 'SULI' 'ECHO' 1
  15.   'TRAC' 'X';
  16.  
  17. GRAPH   = FAUX ;
  18. * GRAPH   = VRAI ;
  19.  
  20. *
  21. *** Methodes possibles :
  22. *
  23. *   'VLH', 'SS'
  24. *
  25.  
  26. METO =  'VLH' ;
  27. * SAFFAC = safety facton on the CFL
  28. SAFFAC = 1.0 ;
  29. * LOGSO = second order reconstruction
  30. LOGSO = VRAI ;
  31. TFINAL = 10. ;
  32. NITER = 40 ;
  33. ZERO = 1.0D-12 ;
  34.  
  35. ************ 
  36. * MAILLAGE *
  37. ************
  38.  
  39. PCEN = 0.5 0.5 ;
  40. A1 = 0.75 0.5 ;
  41. A2 = 1.0 0.5 ;
  42. A3 = 1.0 0.75 ;
  43. A4 = 1.0 1.0 ;
  44. A5 = 0.75 1.0 ;
  45. A6 = 0.5 1.0 ;
  46. A7 = 0.5 0.75 ;
  47.  
  48. MAIL1 = MANU 'QUA4' A1 A3 A5 A7 ;
  49. MAIL1 = MAIL1 ET (MANU 'TRI3' PCEN A1 A7 );
  50. MAIL1 = MAIL1 ET (MANU 'TRI3' A1 A2 A3 );
  51. MAIL1 = MAIL1 ET (MANU 'TRI3' A3 A4 A5 );
  52. MAIL1 = MAIL1 ET (MANU 'TRI3' A5 A6 A7 );
  53.  
  54. MAIL2 = MAIL1;
  55.  
  56. REPE BL1 3 ;
  57.     MAIL2 = MAIL2 TOUR 90 PCEN ;
  58.     MAIL1 = MAIL1 ET MAIL2 ;
  59. FIN BL1 ;
  60.  
  61. ELIM MAIL1 0.001 ;
  62.  
  63. VECX = 1.0 0.0 ;
  64. VECY = 0.0 1.0 ;
  65.  
  66. NELEM = 1 ;
  67. NELEMTOT = (2 * NELEM) + 1 ;
  68.  
  69. MAILTOT = MAIL1 'COUL' 'BLEU' ;
  70.  
  71. IX = (-1 * NELEM) - 1 ;
  72. 'REPE' BL1 NELEMTOT ;
  73.       IX = IX + 1;
  74.       IY =  (-1 * NELEM) - 1 ;
  75.      'REPE' BL2 NELEMTOT ;
  76.             IY = IY + 1; 
  77.             DVECX = IX '*' VECX ;
  78.             DVECY = IY '*' VECY ;
  79.             DVEC = 'PLUS' DVECX DVECY ;
  80.             SI ((IX NEG 0) OU (IY NEG 0)) ;
  81.                MAIL2 = MAIL1 PLUS DVEC ;
  82.                MAILTOT = MAILTOT 'ET' MAIL2 ;
  83.             FINSI ;
  84.       'FIN' BL2 ;
  85. 'FIN' BL1 ;
  86.  
  87. 'ELIM' MAILTOT 0.01 ;
  88. 'ELIM' 0.01 MAIL1 MAILTOT ;
  89.  
  90.  MDNS = 'EULER' ;
  91. $DOMTOT = 'MODE' MAILTOT MDNS  ;
  92. $DOM1   = 'MODE' MAIL1   MDNS  ;
  93. TAB1    = 'DOMA' $DOMTOT 'VF' ;
  94. TAB2    = 'DOMA' $DOM1   'VF' ;
  95. MDOMTOT = TAB1 .  'QUAF' ;
  96. MDOM1   = TAB2 .  'QUAF' ;
  97. 'ELIM' (MDOMTOT 'ET' MDOM1) 1.E-5 ;
  98.  
  99. 'SI' GRAPH ;
  100.  'TRACER' (TAB1 . 'MAILLAGE') 'TITRE' ('CHAINE' 'Maillage ');
  101. 'FINSI' ;   
  102.  
  103.  
  104. ***********************
  105. **** LA TABLE PGAZ ****
  106. ***********************
  107.  
  108. PGAZ = 'TABLE' ;
  109.  
  110. *
  111. **** Ordre des polynoms
  112. *
  113.  
  114. PGAZ . 'NORD' = 4 ;
  115.  
  116. *
  117. **** Especes qui sont dans les equations d'Euler
  118. *
  119.  
  120. PGAZ . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;
  121.  
  122. *
  123. **** Espece qui n'y est pas
  124. *
  125.  
  126.  
  127. PGAZ . 'ESPNEULE' = 'N2  ';
  128.  
  129. *
  130.  
  131. PGAZ .  'H2  ' = 'TABLE'  ;
  132. PGAZ .  'H2O ' = 'TABLE'  ;
  133. PGAZ .  'N2  ' = 'TABLE'  ;
  134. PGAZ .  'O2  ' = 'TABLE'  ;
  135.  
  136. *
  137. **** R (J/Kg/K)
  138. *
  139.  
  140. PGAZ .  'H2  ' . 'R' = 4130.0 ;
  141. PGAZ .  'H2O ' . 'R' = 461.4 ;
  142. PGAZ .  'N2  ' . 'R' = 296.8 ;
  143. PGAZ .  'O2  ' . 'R' = 259.8 ;
  144.  
  145.  
  146. *
  147. **** Regressions polynomials
  148. *
  149.  
  150. PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
  151.                                -2.37281455E-07 1.84701105E-11 ;
  152. PGAZ .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
  153.                               -1.82753232E-08 2.44485692E-12 ;
  154. PGAZ .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
  155.                               8.78233606E-09 -3.05514485E-13 ;
  156. PGAZ .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865
  157.                               2.33636971E-08 -1.53304905E-12;
  158.  
  159. *
  160. **** "Enthalpies" (ou energies) de formations a OK (J/Kg)
  161. *
  162.  
  163. PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
  164. PGAZ .  'H2O ' . 'H0K' = -1.395D7 ;
  165. PGAZ .  'N2  ' . 'H0K' = -2.953D5 ;
  166. PGAZ .  'O2  ' . 'H0K' = -2.634D5 ;
  167.  
  168. *
  169. *** Fin PGAZ
  170. *
  171. *
  172. **** ETAT CENTRE
  173. *
  174.  
  175. rog   = 1.0  ;
  176. ung   = 0.0 ;
  177. utg   = 0.0 ;
  178. retg  = .4291145555695540D+04 ;
  179. yh2g = 0.01 ;
  180. yo2g = 0.2 ;
  181. yh2og = 0.3 ;
  182.  
  183.  
  184. rouxg  = ung '*' rog ;
  185. rouyg  = utg '*' rog ;
  186.  
  187.  
  188. ROC1 = KCHT $DOM1 'SCAL' 'CENTRE' rog ;
  189. ROVC1 = KCHT $DOM1 'VECT' 'CENTRE' (rouxg rouyg );
  190. RETC1 = KCHT $DOM1 'SCAL' 'CENTRE' retg ;
  191. RYH2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2  ' (rog*yh2g) ;
  192. RYO2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'O2  ' (rog*yo2g) ;
  193. RYH2OC1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rog*yh2og) ;
  194.  
  195.  
  196. *
  197. *** Le reste
  198. *
  199.  
  200. rod   =  1.250D-1 ;
  201. und   =  0.0D0 ;
  202. utd   =  0.0D0 ;
  203. retd  = .3598345082089522D+01 ;
  204. yh2d = 0.2 ;
  205. yo2d = 0.4 ;
  206. yh2od = 0.1 ;
  207.  
  208. rouxd  = und '*' rod ;
  209. rouyd  = utd '*' rod ;
  210.  
  211. ROC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' rod ;
  212. ROVC2 = KCHT $DOMTOT 'VECT' 'CENTRE' (rouxd rouyd );
  213. RETC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' retd ;
  214. RYH2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2  '  (rod * yh2d) ;
  215. RYO2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2  '  (rod * yo2d) ;
  216. RYH2OC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rod * yh2od) ;
  217.  
  218. RN  = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' ROC2 ROC1 ;
  219. GN  = 'KCHT'  $DOMTOT 'VECT' 'CENTRE' ROVC2 ROVC1 ;
  220. RETN = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' RETC2 RETC1 ;
  221.  
  222. RYH2  =  KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2  ' RYH2C2 RYH2C1 ;
  223. RYO2  =  KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2  ' RYO2C2 RYO2C1 ;
  224. RYH2O =  KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' RYH2OC2 RYH2OC1 ;
  225.  
  226. RYN =  RYH2 'ET' RYO2 'ET' RYH2O ;
  227.  
  228. ********************************************************
  229. *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM   ***
  230. ********************************************************
  231.  
  232. MOD1     =  'MODELISER' (TAB1 . 'MAILLAGE') 'THERMIQUE' ;
  233.  
  234.  
  235. *
  236. *** GRAPHIQUE DES C.I.
  237. *
  238.  
  239. 'SI' GRAPH ;
  240. * 'SI' FAUX ;
  241. *
  242. *** CREATION DE CHAMELEM
  243. *
  244.  
  245.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  246.   CHM_GN  =  'KCHA' $DOMTOT 'CHAM' GN ;
  247.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  248.   CHM_RYN   =  'KCHA' $DOMTOT 'CHAM' RYN ;
  249.   TRAC CHM_RN MOD1 'TITR'  ('CHAINE'  'RN at t=' 0.0);
  250.   TRAC CHM_GN MOD1 'TITR'  ('CHAINE'  'GN at t=' 0.0);
  251.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' 0.0);
  252.   TRAC CHM_RYN MOD1 'TITR'  ('CHAINE'  'RYN  at t=' 0.0);
  253.  
  254.  
  255. 'FINSI' ;
  256.  
  257. **********************************************************
  258.  
  259.  
  260.  RN0 = 'COPIER' RN ;
  261.  GN0 = 'COPIER' GN ;
  262.  RETN0 = 'COPIER' RETN ;
  263.  RYN0 = 'COPIER' RYN ;
  264.  
  265. *
  266. **** Parameter for the time loop 
  267. *
  268.  
  269. * Names of the residuum components
  270.  
  271.  LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ;
  272.  NOMRN = 'MOTS' 'SCAL' ;
  273.  NOMVN = 'MOTS' 'UX' 'UY';
  274.  NOMY  = 'MOTS' 'H2' 'O2' 'H2O' ;
  275.  
  276. *
  277. **** Geometrical coefficient to compute gradients
  278. *
  279.  
  280.  GRADR CACCA COEFSCAL =  'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  281.      NOMRN     RN ;
  282.  
  283.  GRADV CACCA COEFVECT =  'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  284.      NOMVN     GN ;
  285.  
  286.  TPS = 0.0 ;
  287.  
  288.  
  289. *
  290. **** Temporal loop
  291. *
  292.  
  293.  'MESSAGE' ;
  294.  'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
  295.  'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
  296.  'MESSAGE' ;
  297.  
  298.   VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
  299.  
  300.   TNM1 = 'COPIER' TN ;
  301.  
  302.  'TEMPS' 'ZERO' ;
  303.  'REPETER' BL1 NITER ;
  304.  
  305. *
  306. **** Primitive variables
  307. *
  308.  
  309.     VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
  310.       TNM1  ;
  311.  
  312.     TNM1 = 'COPIER' TN ;
  313.  
  314.     'SI' LOGSO ;
  315.  
  316.        GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  317.          NOMRN  RN  'GRADGEO' COEFSCAL ;
  318.  
  319.        GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  320.          NOMRN  PN  'GRADGEO' COEFSCAL ;
  321.  
  322.        GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  323.          NOMVN  VN   'GRADGEO'  COEFVECT ;
  324.  
  325.        GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  326.          NOMY  YN  'GRADGEO' COEFSCAL ;
  327.  
  328.        ROF VITF PF YF    = 'PRET' 'PERFTEMP'  2 1
  329.                             $DOMTOT  PGAZ
  330.                             RN    GRADR  ALR
  331.                             VN    GRADV  ALV
  332.                             PN    GRADP  ALP
  333.                             YN    GRADY  ALY ;
  334.  
  335.     'SINON' ;
  336.  
  337.        ROF VITF PF YF    = 'PRET' 'PERFTEMP'  1 1 
  338.                               $DOMTOT  PGAZ
  339.                               RN   
  340.                               VN   
  341.                               PN   
  342.                               YN    ;
  343.  
  344.     'FINSI' ;
  345.  
  346.     RESIDU DELTAT =  'KONV' 'VF' 'PERFTEMP' 'RESI' METO
  347.          $DOMTOT PGAZ LISTINCO  ROF VITF PF YF ;
  348.  
  349.     DT_CON = SAFFAC '*' DELTAT ;     
  350. *
  351. **** The time step linked to tps
  352. *
  353.  
  354.     DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
  355.  
  356. *
  357. **** Total time step
  358. *
  359.  
  360.     DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
  361.  
  362. *
  363. **** Increment of the variables (convection)
  364. *
  365.  
  366.     RESIDU = DTMIN '*' RESIDU ;
  367.  
  368.     DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
  369.     DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
  370.     DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
  371.     DRYN  = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
  372.             ('MOTS' 'H2' 'O2' 'H2O') ;
  373.  
  374.     TPS  = TPS '+' DTMIN ;
  375.     RN   = RN '+' DRN ;
  376.     GN   = GN '+' DGN ;
  377.     RETN = RETN '+' DRETN ;
  378.     RYN = RYN '+' DRYN ;
  379.  
  380.     'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
  381.         'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
  382.     'FINSI' ;
  383.  
  384.     'SI' (TPS > TFINAL) ;
  385.        'QUITTER' BL1 ;
  386.     'FINSI' ;
  387.  
  388.  'FIN' BL1 ;
  389.  'TEMPS' ;
  390.  
  391. **********************************************************
  392.  
  393. *
  394. **** Les variables primitives
  395. *
  396.  
  397. VN PN TN YN GAMN =   'PRIM' 'PERFTEMP' PGAZ
  398.            RN GN RETN RYN ;
  399.  
  400. CN = (GAMN * (PN / RN)) '**' 0.5 ;
  401.  
  402. *
  403. *** GRAPHIQUE DES SOLUTIONS
  404. *
  405.  
  406. 'SI' GRAPH ;
  407.  
  408. *
  409. *** CREATION DE CHAMELEM
  410. *
  411.  
  412.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  413.   CHM_GN   =  'KCHA' $DOMTOT 'CHAM' GN ;
  414.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  415.   CHM_RYN  =  'KCHA' $DOMTOT 'CHAM' RYN ;  
  416.   TRAC CHM_RN  MOD1 'TITR'  ('CHAINE'  'RO at t=' TFIN);
  417.   TRAC CHM_GN  MOD1 'TITR'  ('CHAINE'  'GN at t=' TFIN);
  418.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' TFIN);
  419.   TRAC CHM_RYN MOD1 'TITR'  ('CHAINE'  'RYN at t=' TFIN);
  420.  
  421. 'FINSI' ;
  422.  
  423. * Tests: invariance par rotation
  424.  
  425. RN1 = COPIER RN ;
  426. GN1 = COPI GN ;
  427. RETN1 = COPIER RETN ;
  428. RYN1 = COPIER RYN ;
  429.  
  430. MAIL0 = EXTR RN 'MAILLAGE' ;
  431.  
  432. REPETER BL1 3 ;
  433.  
  434.    RN1 = RN1 TOUR 90 PCEN ;
  435.    GN1 = GN1 TOUR 90 PCEN ;
  436.    RETN1 = RETN1 TOUR 90  PCEN ;
  437.    RYN1 = RYN1 TOUR 90 PCEN ;
  438.  
  439.    MAIL1 = EXTR RN1 'MAILLAGE' ;
  440.    ELIM 0.001 MAIL0 MAIL1 ;
  441.    ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ;
  442.    SI (ERR1 > 1.0D-10);
  443.       ERRE 5;
  444.    FINSI ;
  445.  
  446.    MAIL1 = EXTR GN1 'MAILLAGE' ;
  447.    ELIM 0.001 MAIL0 MAIL1 ;
  448.    ERR1 = MAXI ((GN1 - GN) / (RN * CN)) 'ABS' ;
  449.    SI (ERR1 > 1.0D-10);
  450.       ERRE 5;
  451.    FINSI ;
  452.  
  453.    MAIL1 = EXTR RETN1 'MAILLAGE' ;
  454.    ELIM 0.001 MAIL0 MAIL1 ;
  455.    ERR1 = MAXI ((RETN1 - RETN) / RETN) 'ABS' ;
  456.    SI (ERR1 > 1.0D-10);
  457.       ERRE 5;
  458.    FINSI ;
  459.  
  460.    MAIL1 = EXTR RYN1 'MAILLAGE' ;
  461.    ELIM 0.001 MAIL0 MAIL1 ;
  462.    ERR1 = MAXI ((RYN1 - RYN) / RN) 'ABS' ;
  463.    SI (ERR1 > 1.0D-10);
  464.       ERRE 5;
  465.    FINSI ;
  466.  
  467. FIN BL1 ;
  468.  
  469. * Symetrie
  470.  
  471. MAIL1 RN1 = MAIL0 RN SYME 'DROI' (0 0.5) (1.0 0.5) ;
  472. ELIM 0.001 MAIL0 MAIL1 ;
  473. ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ;
  474. SI (ERR1 > 1.0D-10);
  475.       ERRE 5;
  476. FINSI ;
  477.  
  478. FIN ;
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  

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