Télécharger la plaquette Cast3M


fr en

Télécharger carre_calper.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  carre_calper.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ******************************************************************
  5. * CALCUL DU QUARRE A CHOC                                        *
  6. * Modele: gaz multi-especes "thermally perfect"                  *
  7. * Gaz multi-especes "calorically 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.  
  100. 'SI' GRAPH ;
  101.  'TRACER' (TAB1 . 'MAILLAGE') 'TITRE' ('CHAINE' 'Maillage ');
  102. 'FINSI' ;   
  103.  
  104. ***********************
  105. **** LA TABLE PGAZ ****
  106. ***********************
  107.  
  108. PGAZ = 'TABLE' ;
  109.  
  110. *
  111. **** Ordre des polynoms cv_i
  112. *
  113.  
  114. PGAZ . 'NORD' = 0 ;
  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.  
  151. *
  152. * Cas particulier: gaz "calorically perfect"
  153. *
  154.  
  155. PGAZ .  'H2  ' . 'A' = 'PROG' .14571861D+05 ; 
  156. PGAZ .  'H2O ' . 'A' = 'PROG' .26589930D+04 ; 
  157. PGAZ .  'N2  ' . 'A' = 'PROG' .10024563D+04 ; 
  158. PGAZ .  'O2  ' . 'A' = 'PROG' .92885670D+03 ; 
  159.  
  160. *
  161. **** "Enthalpies" (ou energies) de formations a OK (J/Kg)
  162. *     Note: ce sont des entites numeriques
  163. *     h_i = h_i(T0) '-' \int_0^{T0} cp_i(x) dx =
  164. *           h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
  165. *
  166. *     Pour H2, H20, O2, N2 on a:
  167. *
  168. *     T0 = 298.15
  169. *     
  170.  
  171. PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
  172. PGAZ .  'H2O ' . 'H0K' = -1.395D7 ;
  173. PGAZ .  'N2  ' . 'H0K' = -2.953D5 ;
  174. PGAZ .  'O2  ' . 'H0K' = -2.634D5 ;
  175.  
  176. *
  177. *** Fin PGAZ
  178. *
  179.  
  180.  
  181. *
  182. **** ETAT CENTRE
  183. *
  184.  
  185. rog   = 1.0  ;
  186. ung   = 0.0 ;
  187. utg   = 0.0 ;
  188. retg  = .4291145555695540D+04 ;
  189. yh2g = 0.01 ;
  190. yo2g = 0.2 ;
  191. yh2og = 0.3 ;
  192.  
  193.  
  194. rouxg  = ung '*' rog ;
  195. rouyg  = utg '*' rog ;
  196.  
  197.  
  198. ROC1 = KCHT $DOM1 'SCAL' 'CENTRE' rog ;
  199. ROVC1 = KCHT $DOM1 'VECT' 'CENTRE' (rouxg rouyg );
  200. RETC1 = KCHT $DOM1 'SCAL' 'CENTRE' retg ;
  201. RYH2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2  ' (rog*yh2g) ;
  202. RYO2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'O2  ' (rog*yo2g) ;
  203. RYH2OC1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rog*yh2og) ;
  204.  
  205.  
  206. *
  207. *** Le reste
  208. *
  209.  
  210. rod   =  1.250D-1 ;
  211. und   =  0.0D0 ;
  212. utd   =  0.0D0 ;
  213. retd  = .3598345082089522D+01 ;
  214. yh2d = 0.2 ;
  215. yo2d = 0.4 ;
  216. yh2od = 0.1 ;
  217.  
  218. rouxd  = und '*' rod ;
  219. rouyd  = utd '*' rod ;
  220.  
  221. ROC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' rod ;
  222. ROVC2 = KCHT $DOMTOT 'VECT' 'CENTRE' (rouxd rouyd );
  223. RETC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' retd ;
  224. RYH2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2  '  (rod * yh2d) ;
  225. RYO2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2  '  (rod * yo2d) ;
  226. RYH2OC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rod * yh2od) ;
  227.  
  228. RN  = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' ROC2 ROC1 ;
  229. GN  = 'KCHT'  $DOMTOT 'VECT' 'CENTRE' ROVC2 ROVC1 ;
  230. RETN = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' RETC2 RETC1 ;
  231.  
  232. RYH2  =  KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2  ' RYH2C2 RYH2C1 ;
  233. RYO2  =  KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2  ' RYO2C2 RYO2C1 ;
  234. RYH2O =  KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' RYH2OC2 RYH2OC1 ;
  235.  
  236. RYN =  RYH2 'ET' RYO2 'ET' RYH2O ;
  237.  
  238. ********************************************************
  239. *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM   ***
  240. ********************************************************
  241.  
  242. MOD1     =  'MODELISER'  (TAB1 . 'MAILLAGE') 'THERMIQUE' ;
  243.  
  244.  
  245. *
  246. *** GRAPHIQUE DES C.I.
  247. *
  248.  
  249. 'SI' GRAPH ;
  250. * 'SI' FAUX ;
  251. *
  252. *** CREATION DE CHAMELEM
  253. *
  254.  
  255.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  256.   CHM_GN  =  'KCHA' $DOMTOT 'CHAM' GN ;
  257.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  258.   CHM_RYN   =  'KCHA' $DOMTOT 'CHAM' RYN ;
  259.   TRAC CHM_RN MOD1 'TITR'  ('CHAINE'  'RN at t=' 0.0);
  260.   TRAC CHM_GN MOD1 'TITR'  ('CHAINE'  'GN at t=' 0.0);
  261.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' 0.0);
  262.   TRAC CHM_RYN MOD1 'TITR'  ('CHAINE'  'RYN  at t=' 0.0);
  263.  
  264.  
  265. 'FINSI' ;
  266.  
  267. **********************************************************
  268.  
  269.  
  270.  RN0 = 'COPIER' RN ;
  271.  GN0 = 'COPIER' GN ;
  272.  RETN0 = 'COPIER' RETN ;
  273.  RYN0 = 'COPIER' RYN ;
  274.  
  275. *
  276. **** Parameter for the time loop 
  277. *
  278.  
  279. * Names of the residuum components
  280.  
  281.  LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ;
  282.  NOMRN = 'MOTS' 'SCAL' ;
  283.  NOMVN = 'MOTS' 'UX' 'UY' ;
  284.  NOMY  = 'MOTS' 'H2' 'O2' 'H2O' ;
  285.  
  286. *
  287. **** Geometrical coefficient to compute gradients
  288. *
  289.  
  290.  GRADR CACCA COEFSCAL =  'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  291.         NOMRN  RN ;
  292.  
  293.  GRADV CACCA COEFVECT =  'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  294.         NOMVN  GN ;
  295.  
  296.  TPS = 0.0 ;
  297.  
  298.  
  299. *
  300. **** Temporal loop
  301. *
  302.  
  303.  'MESSAGE' ;
  304.  'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
  305.  'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
  306.  'MESSAGE' ;
  307.  
  308.   VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
  309.  
  310.   TNM1 = 'COPIER' TN ;
  311.  
  312.  'TEMPS' 'ZERO' ;
  313.  'REPETER' BL1 NITER ;
  314.  
  315. *
  316. **** Primitive variables
  317. *
  318.  
  319.     VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
  320.       TNM1  ;
  321.  
  322.     TNM1 = 'COPIER' TN ;
  323.  
  324.     'SI' LOGSO ;
  325.  
  326.        GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  327.          NOMRN RN  'GRADGEO' COEFSCAL ;
  328.  
  329.        GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  330.          NOMRN PN  'GRADGEO' COEFSCAL ;
  331.  
  332.        GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  333.          NOMVN  VN   'GRADGEO'  COEFVECT ;
  334.  
  335.        GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  336.          NOMY  YN  'GRADGEO' COEFSCAL ;
  337.  
  338.        ROF VITF PF YF    = 'PRET' 'PERFTEMP'  2 1
  339.                             $DOMTOT  PGAZ
  340.                             RN    GRADR  ALR
  341.                             VN    GRADV  ALV
  342.                             PN    GRADP  ALP
  343.                             YN    GRADY  ALY ;
  344.  
  345.     'SINON' ;
  346.  
  347.        ROF VITF PF YF    = 'PRET' 'PERFTEMP'  1 1 
  348.                               $DOMTOT  PGAZ
  349.                               RN   
  350.                               VN   
  351.                               PN   
  352.                               YN    ;
  353.  
  354.     'FINSI' ;
  355.  
  356.     RESIDU DELTAT =  'KONV' 'VF' 'PERFTEMP' 'RESI' METO
  357.          $DOMTOT PGAZ LISTINCO  ROF VITF PF YF ;
  358.  
  359.     DT_CON = SAFFAC '*' DELTAT ;     
  360. *
  361. **** The time step linked to tps
  362. *
  363.  
  364.     DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
  365.  
  366. *
  367. **** Total time step
  368. *
  369.  
  370.     DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
  371.  
  372. *
  373. **** Increment of the variables (convection)
  374. *
  375.  
  376.     RESIDU = DTMIN '*' RESIDU ;
  377.  
  378.     DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
  379.     DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
  380.     DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
  381.     DRYN  = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
  382.             ('MOTS' 'H2' 'O2' 'H2O') ;
  383.  
  384.     TPS  = TPS '+' DTMIN ;
  385.     RN   = RN '+' DRN ;
  386.     GN   = GN '+' DGN ;
  387.     RETN = RETN '+' DRETN ;
  388.     RYN = RYN '+' DRYN ;
  389.  
  390.     'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
  391.         'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
  392.     'FINSI' ;
  393.  
  394.     'SI' (TPS > TFINAL) ;
  395.        'QUITTER' BL1 ;
  396.     'FINSI' ;
  397.  
  398.  'FIN' BL1 ;
  399.  'TEMPS' ;
  400.  
  401. **********************************************************
  402.  
  403. *
  404. **** Les variables primitives
  405. *
  406.  
  407. VN PN TN YN GAMN =   'PRIM' 'PERFTEMP' PGAZ
  408.            RN GN RETN RYN ;
  409.  
  410. CN = (GAMN * (PN / RN)) '**' 0.5 ;
  411.  
  412. *
  413. *** GRAPHIQUE DES SOLUTIONS
  414. *
  415.  
  416. 'SI' GRAPH ;
  417.  
  418. *
  419. *** CREATION DE CHAMELEM
  420. *
  421.  
  422.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  423.   CHM_GN   =  'KCHA' $DOMTOT 'CHAM' GN ;
  424.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  425.   CHM_RYN  =  'KCHA' $DOMTOT 'CHAM' RYN ;  
  426.   TRAC CHM_RN  MOD1 'TITR'  ('CHAINE'  'RO at t=' TFIN);
  427.   TRAC CHM_GN  MOD1 'TITR'  ('CHAINE'  'GN at t=' TFIN);
  428.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' TFIN);
  429.   TRAC CHM_RYN MOD1 'TITR'  ('CHAINE'  'RYN at t=' TFIN);
  430.  
  431. 'FINSI' ;
  432.  
  433. * Tests: invariance par rotation
  434.  
  435. RN1 = COPIER RN ;
  436. GN1 = COPI GN ;
  437. RETN1 = COPIER RETN ;
  438. RYN1 = COPIER RYN ;
  439.  
  440. MAIL0 = EXTR RN 'MAILLAGE' ;
  441.  
  442. REPETER BL1 3 ;
  443.  
  444.    RN1 = RN1 TOUR 90 PCEN ;
  445.    GN1 = GN1 TOUR 90 PCEN ;
  446.    RETN1 = RETN1 TOUR 90  PCEN ;
  447.    RYN1 = RYN1 TOUR 90 PCEN ;
  448.  
  449.    MAIL1 = EXTR RN1 'MAILLAGE' ;
  450.    ELIM 0.001 MAIL0 MAIL1 ;
  451.    ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ;
  452.    SI (ERR1 > 1.0D-10);
  453.       ERRE 5;
  454.    FINSI ;
  455.  
  456.    MAIL1 = EXTR GN1 'MAILLAGE' ;
  457.    ELIM 0.001 MAIL0 MAIL1 ;
  458.    ERR1 = MAXI ((GN1 - GN) / (RN * CN)) 'ABS' ;
  459.    SI (ERR1 > 1.0D-10);
  460.       ERRE 5;
  461.    FINSI ;
  462.  
  463.    MAIL1 = EXTR RETN1 'MAILLAGE' ;
  464.    ELIM 0.001 MAIL0 MAIL1 ;
  465.    ERR1 = MAXI ((RETN1 - RETN) / RETN) 'ABS' ;
  466.    SI (ERR1 > 1.0D-10);
  467.       ERRE 5;
  468.    FINSI ;
  469.  
  470.    MAIL1 = EXTR RYN1 'MAILLAGE' ;
  471.    ELIM 0.001 MAIL0 MAIL1 ;
  472.    ERR1 = MAXI ((RYN1 - RYN) / RN) 'ABS' ;
  473.    SI (ERR1 > 1.0D-10);
  474.       ERRE 5;
  475.    FINSI ;
  476.  
  477. FIN BL1 ;
  478.  
  479. * Symetrie
  480.  
  481. MAIL1 RN1 = MAIL0 RN SYME 'DROI' (0 0.5) (1.0 0.5) ;
  482. ELIM 0.001 MAIL0 MAIL1 ;
  483. ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ;
  484. SI (ERR1 > 1.0D-10);
  485.       ERRE 5;
  486. FINSI ;
  487.  
  488. FIN ;
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  
  507.  
  508.  
  509.  

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