Télécharger la plaquette Cast3M


fr en

Télécharger prim_ther_mono_3D.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  prim_ther_mono_3D.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***********************************************************
  5. **** APPROCHE VF "Cell-Centred Formulation" pour la    ****
  6. **** solution des                                      **** 
  7. **** Equations d'Euler pour un gaz parfait             ****
  8. **** OPERATEUR  PRIM                                   ****
  9. **** Gaz monoespece "thermally perfect"                ****
  10. **** Cas 3D                                            ****
  11. ****                                                   ****
  12. **** A. BECCANTINI DRN/DMT/SEMT/TTMF  JUILLET 1998     ****
  13. ***********************************************************
  14. *
  15. 'OPTION'  'DIME' 3 ;
  16. 'OPTI' 'ELEM' 'CUB8' ;
  17. 'OPTION'  'ECHO'  0 ;
  18. 'OPTION'  'TRAC' 'X';
  19.  
  20. *
  21. *
  22. **** GRAPH
  23. *
  24. *
  25. GRAPH = FAUX ;
  26. * GRAPH = VRAI ;
  27. *
  28.  
  29. ERRMAX = 1.0D-6 ;        
  30.  
  31. 'MESSAGE' ;
  32. 'MESSAGE' ('cv = polynom de 4-eme degree');
  33. 'MESSAGE' ;
  34.  
  35. *
  36. *** On considere une melange H2, O2, H2O, N2
  37. *
  38. * De la table JANAF on a, pour les cp (J/Kg/K) 
  39. *
  40.  
  41. LTEMP = 'PROG' .2000D+03   .3000D+03   .4000D+03   .5000D+03   .6000D+03   
  42.          .7000D+03   .8000D+03   .9000D+03   .1000D+04   .1100D+04   
  43.          .1200D+04   .1300D+04   .1400D+04   .1500D+04   .1600D+04   
  44.          .1700D+04   .1800D+04   .1900D+04   .2000D+04   .2100D+04   
  45.          .2200D+04   .2300D+04   .2400D+04   .2500D+04   .2600D+04 ;
  46.  
  47. LTEMP = LTEMP 'ET' ('PROG'                     
  48.          .2700D+04   .2800D+04   .2900D+04   .3000D+04   .3100D+04
  49.          .3200D+04   .3300D+04   .3400D+04   .3500D+04   .3600D+04   
  50.          .4000D+04   .4500D+04   .5000D+04   .5500D+04   .6000D+04) ;
  51.  
  52. LCVH2 = 'PROG' .9944D+04   .1018D+05   .1035D+05   .1039D+05   .1042D+05   
  53.          .1048D+05   .1057D+05   .1070D+05   .1086D+05   .1104D+05   
  54.          .1125D+05   .1147D+05   .1168D+05   .1190D+05   .1211D+05   
  55.          .1231D+05   .1251D+05   .1270D+05   .1288D+05   .1305D+05   
  56.          .1321D+05   .1337D+05   .1351D+05   .1365D+05   .1378D+05 ;
  57.  
  58. LCVH2 = LCVH2 'ET' ('PROG'                   
  59.          .1391D+05   .1403D+05   .1415D+05   .1426D+05   .1438D+05 
  60.          .1449D+05   .1460D+05   .1471D+05   .1480D+05   .1490D+05
  61.          .1528D+05   .1573D+05   .1613D+05   .1646D+05   .1669D+05 ) ;
  62.  
  63.  
  64.  
  65. ***********************
  66. **** LA TABLE PGAZ ****
  67. ***********************
  68.  
  69. PGAZ = 'TABLE' ;
  70.  
  71.  
  72. PGAZ . 'NORD' = 4 ;
  73. *
  74. **** La seul espece 
  75. *
  76.  
  77. PGAZ . 'ESPNEULE' = 'H2  ';
  78.  
  79. *
  80.  
  81. PGAZ .  'H2  ' = 'TABLE'  ;
  82. *
  83. **** R (J/Kg/K)
  84. *
  85.  
  86. PGAZ .  'H2  ' . 'R' = 4130.0 ;
  87.  
  88. *
  89. **** Regressions polynomials
  90. *
  91.  
  92.    PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
  93.                                -2.37281455E-07 1.84701105E-11 ;
  94.  
  95. *
  96. **** "Enthalpies" (ou energies) de formations a OK (J/Kg)
  97. *     Note: ce sont des entites numeriques
  98.  
  99. PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
  100.  
  101. *
  102. *** Fin PGAZ
  103. *
  104.  
  105. A0H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 1 ;
  106. A1H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 2 ;
  107. A2H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 3 ;
  108. A3H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 4 ;
  109. A4H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 5 ;
  110.  
  111.   LCV1H2 = 'PROG' ;                                                            
  112.   LCV1H2O = 'PROG' ;                                                            
  113.   LCV1N2 = 'PROG' ;                                                            
  114.     LCV1O2 = 'PROG' ;                                                            
  115.  
  116.   'REPETER' BL1 ('DIME' LTEMP ) ;
  117.      T = 'EXTRAIRE' LTEMP &BL1 ;
  118.      T2 = T * T ;
  119.      T3 = T2 * T ;
  120.      T4 = T3 * T;
  121.      LCV1H2 = LCV1H2 'ET' ('PROG'
  122.         (A0H2 '+' (A1H2 * T) '+' (A2H2 * T2) '+' (A3H2 * T3) '+'
  123.         (A4H2 * T4))) ;
  124.   'FIN' BL1 ;
  125.  
  126.   EVCVH2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCVH2 ;
  127.   EVCV1H2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCV1H2 ;
  128.  
  129.   TAB1 = 'TABLE'  ;                                                             
  130.   TAB1 . 'TITRE' = 'TABLE'  ;                                                   
  131.   TAB1 . 1 ='TIRR ';                                                         
  132.   TAB1 . 'TITRE' . 1  =  'JANAF DATA';                                  
  133.   TAB1 . 2 = 'MARQ CROI';                                                     
  134.   TAB1 . 'TITRE' . 2  =  'POLYN. REGR.';                                  
  135.  
  136.   'SI' GRAPH ;
  137.      'DESSIN'  (EVCVH2 'ET' EVCV1H2) 'LEGE'  TAB1 'TITRE' 'H2' ;
  138.   'FINSI' ;
  139.  
  140.  
  141. *************************** 
  142. ***** DOMAINE SPATIAL  ****
  143. ***************************
  144.  
  145. ************ 
  146. * MAILLAGE *
  147. ************
  148.  
  149. NX =  10 ;
  150. NY =  2  ;
  151. NZ =  2 ;
  152.  
  153. L  = 1.0D0  ;
  154. DX = L '/' NX '/' 2.0D0 ;
  155. H  = NY '*' DX  ;
  156. P  = NZ '*' DX ;
  157. xD =  0.5D0 '*' L ;
  158. xG = -1.0D0 '*' xD  ;
  159. yH =  0.5D0 '*' H   ;
  160. yB =  -1.0D0 '*' yH ;
  161. zV = 0.5D0 '*' P    ;
  162. zR = -1.0D0 '*' zV  ;
  163.  
  164. A1 = xG yB zR ;
  165. A2 = 0.0D0 yB zR ;
  166. A7 = xG yB  zV ;
  167. A8 = 0.0D0 yB  zV ;
  168.  
  169.  
  170. BAS1  = A1 'DROI' NX A2 ;
  171. BAS3  = A7 'DROI' NX A8 ;
  172. BAS5  = A1 'DROI' NZ A7 ;
  173. BAS6  = A2 'DROI' NZ A8 ;
  174.  
  175.  
  176. S11 = 'DALL' BAS3 BAS6 BAS1 BAS5 ;
  177.  
  178.  
  179. DOM1   = S11 'VOLU' NZ 'TRANS' (0.0 H 0.0) ;
  180. $DOM1 = 'DOMA' DOM1 ;
  181.  
  182. *
  183. **** Pression, densite, temperature et vitesse
  184. *
  185.  
  186. MAIL =   $DOM1 . 'CENTRE' ;
  187.  
  188. RN   = 'BRUI' 'BLAN' 'UNIF' 1.0  0.1  MAIL ;
  189. TN   = 'BRUI' 'BLAN' 'UNIF' 3900  2000  MAIL ;
  190.  
  191. VNX = 'BRUI' 'BLAN' 'UNIF' 1.0  0.1  MAIL ;
  192. VNY = 'BRUI' 'BLAN' 'UNIF' 2.0  0.1  MAIL ;
  193. VNZ = 'BRUI' 'BLAN' 'UNIF' 3.0  0.1  MAIL ;
  194.  
  195. VN = ('NOMC' VNX 'UX' 'NATU' 'DISCRET' ) 'ET'
  196.      ('NOMC' VNY 'UY' 'NATU' 'DISCRET' ) 'ET'
  197.      ('NOMC' VNZ 'UZ' 'NATU' 'DISCRET' ) ;
  198.  
  199. GN = ('NOMC' (VNX * RN) 'UX' 'NATU' 'DISCRET' ) 'ET'
  200.      ('NOMC' (VNY * RN) 'UY' 'NATU' 'DISCRET' ) 'ET'
  201.      ('NOMC' (VNZ * RN) 'UZ' 'NATU' 'DISCRET' ) ;
  202.  
  203. *
  204. * Calcul de l'energie thermique et du gamma
  205. *
  206. * Ether = \sum_{i=1,4} Y_i \int_{0}^{T} cv_i(x) dx
  207. * gamma = (\sum_{i=1,4} Y_i cp_i) / (\sum_{i=1,4} Y_i cv_i)
  208. *
  209.  
  210. T2 = TN * TN ;
  211. T3 = T2 * TN ;
  212. T4 = T3 * TN ;
  213. T5 = T4 * TN ;
  214.  
  215. ETHER =  ((A0H2 * TN) '+' ((A1H2 '/' 2.0) * T2) '+'
  216.          ((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
  217.          ((A4H2 '/' 5.0) * T5)) ;
  218.  
  219. LISCOM2 = 'MOTS' 'UX  ' 'UY  ' 'UZ  ';
  220. RECIN = 'PSCAL' (0.5 '*' VN) GN LISCOM2 LISCOM2 ;
  221.  
  222. RET = RECIN '+' (ETHER * RN) ;
  223.  
  224.  CVTOT =  (A0H2 '+' (A1H2 * TN) '+'
  225.          (A2H2 * T2) '+' (A3H2 * T3) '+'
  226.          (A4H2 * T4)) ; 
  227.  
  228.  
  229.  CPTOT  = CVTOT  '+' (PGAZ . 'H2  ' . 'R') ;
  230.  
  231.  GAMN = CPTOT '/' CVTOT ;
  232.  
  233.  RTOT = CPTOT '-' CVTOT ;
  234.  
  235.  PN = (RTOT '*' TN) '*' RN ;
  236.  
  237. *
  238. **** PRIM
  239. *
  240.  
  241. VITESSE PRES TEMPNEW GAMNEW
  242.        = 'PRIM' 'PERFTEMP' PGAZ RN GN RET ;
  243.  
  244.  
  245. *
  246. **** TEST
  247. *
  248.  
  249. ERRV = ('MAXIMUM' (VN '-' VITESSE) 'ABS')
  250.         '/' ('MAXIMUM' VN) ;
  251.  
  252. ERRP  = 'MAXIMUM' ((PRES '-' PN)  '/'  PN) 'ABS';
  253.  
  254. ERRT  = 'MAXIMUM' ((TEMPNEW '-' TN) '/' TN) 'ABS';  
  255.  
  256. ERRG = 'MAXIMUM' ((GAMN '-' GAMNEW) '/' GAMN) 'ABS'; 
  257.  
  258.  
  259. 'SI' (('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG  ) '>'
  260.           ERRMAX));
  261.      'MESSAGE' ('CHAINE' 'Erreur maximum');
  262.      'LISTE' ('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG
  263.       ));
  264.      'ERREUR'  5;
  265. 'FINSI' ;
  266.  
  267. VITESSE PRES TEMPNEW2 GAMNEW
  268.        = 'PRIM' 'PERFTEMP' PGAZ RN GN RET TEMPNEW;
  269.  
  270. ERRT  = 'MAXIMUM' ((TEMPNEW '-' TEMPNEW2) '/' TEMPNEW) 'ABS';
  271.  
  272. 'SI' (ERRT > ERRMAX);
  273.      'MESSAGE' ('CHAINE' 'Erreur maximum' ERRT) ;
  274.      'ERREUR'  5;
  275. 'FINSI' ;
  276.  
  277. **** On teste aussi la difference entre l'energie initiale et celle la.
  278.  
  279. T2 = TEMPNEW * TEMPNEW ;
  280. T3 = T2 * TEMPNEW ;
  281. T4 = T3 * TEMPNEW ;
  282. T5 = T4 * TEMPNEW ;
  283.  
  284. ETHER_N =  ((A0H2 * TEMPNEW) '+' ((A1H2 '/' 2.0) * T2) '+'
  285.          ((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
  286.          ((A4H2 '/' 5.0) * T5)) ;
  287.  
  288. ERRET = 'MAXIMUM' ((ETHER '-' ETHER_N) '/' ETHER) 'ABS' ;
  289.  
  290.  
  291. 'SI' (ERRET > ERRMAX);
  292.      'MESSAGE' ('CHAINE' 'Erreur maximum' ERRET) ;
  293.      'ERREUR'  5;
  294. 'FINSI' ;
  295.  
  296.  
  297. *) Degre de polynome = 0
  298.  
  299. ERRMAX = 1.0D-10 ;         
  300.  
  301. 'MESSAGE' ;
  302. 'MESSAGE' ('cv = polynom constant');
  303. 'MESSAGE' ;
  304.  
  305.  
  306. ***********************
  307. **** LA TABLE PGAZ ****
  308. ***********************
  309.  
  310. PGAZ = 'TABLE' ;
  311.  
  312. *
  313. **** Degre des polynoms
  314. *
  315.  
  316. PGAZ . 'NORD' = 0 ;
  317.  
  318. *
  319. **** Especes qui sont dans les equations d'Euler
  320. *
  321.  
  322. *
  323. **** Espece qui n'y est pas
  324. *
  325.  
  326.  
  327. PGAZ . 'ESPNEULE' = 'H2  ';
  328.  
  329. *
  330.  
  331. PGAZ .  'H2  ' = 'TABLE'  ;
  332.  
  333. *
  334. **** R (J/Kg/K)
  335. *
  336.  
  337. PGAZ .  'H2  ' . 'R' = 4130.0 ;
  338.  
  339.  
  340. *
  341. **** Regressions polynomials
  342. *
  343.  
  344.  
  345. *
  346. * Cas particulier: gaz "calorically perfect"
  347. *
  348.  
  349. PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 ;
  350.  
  351. *
  352. **** "Enthalpies" (ou energies) de formations a OK (J/Kg)
  353. *
  354.  
  355. PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
  356.  
  357. *
  358. *** Fin PGAZ
  359. *
  360.  
  361. A0H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 1 ;
  362.  
  363.  
  364. *
  365. **** Pression, densite, temperature et vitesse
  366. *
  367.  
  368. RN   = 'BRUI' 'BLAN' 'UNIF' 1.0  0.1  MAIL ;
  369. TN   = 'BRUI' 'BLAN' 'UNIF' 3900  200  MAIL ;
  370. YH2  = 'BRUI' 'BLAN' 'UNIF' 0.3  0.1  MAIL ;
  371. YO2  = 'BRUI' 'BLAN' 'UNIF' 0.2  0.1  MAIL ;
  372. YH2O  = 'BRUI' 'BLAN' 'UNIF' 0.1  0.05  MAIL ;
  373.  
  374. YN2 = 1.0 '-' (YH2 '+' YO2 '+' YH2O) ;
  375.  
  376. RYH2 = RN * YH2 ;
  377. RYH2O = RN * YH2O ;
  378. RYO2 = RN * YO2 ;
  379. RYN2 = RN * YN2 ;
  380.  
  381. YH2  = 'NOMC' YH2  'H2  ' 'NATU' 'DISCRET' ;
  382. YH2O = 'NOMC' YH2O 'H2O ' 'NATU' 'DISCRET' ;
  383. YO2  = 'NOMC' YO2  'O2  ' 'NATU' 'DISCRET' ;
  384. YN2 = 'NOMC' YN2 'N2  ' 'NATU' 'DISCRET' ;
  385. YN = YH2 'ET' YO2 'ET' YH2O ; 
  386.  
  387. RYH2 =  'NOMC' RYH2  'H2  ' 'NATU' 'DISCRET' ;
  388. RYH2O = 'NOMC' RYH2O 'H2O ' 'NATU' 'DISCRET' ;
  389. RYO2 =  'NOMC' RYO2  'O2  ' 'NATU' 'DISCRET' ;
  390. RYN2 =  'NOMC' RYN2  'N2  ' 'NATU' 'DISCRET' ;
  391.  
  392. VNX = 'BRUI' 'BLAN' 'UNIF' 1.0  0.1  MAIL ;
  393. VNY = 'BRUI' 'BLAN' 'UNIF' 2.0  0.1  MAIL ;
  394. VNZ = 'BRUI' 'BLAN' 'UNIF' 3.0  0.1  MAIL ;
  395.  
  396. VN = ('NOMC' VNX 'UX' 'NATU' 'DISCRET' ) 'ET'
  397.      ('NOMC' VNY 'UY' 'NATU' 'DISCRET' ) 'ET'
  398.      ('NOMC' VNZ 'UZ' 'NATU' 'DISCRET' ) ;
  399.  
  400. GN = ('NOMC' (VNX * RN) 'UX' 'NATU' 'DISCRET' ) 'ET'
  401.      ('NOMC' (VNY * RN) 'UY' 'NATU' 'DISCRET' ) 'ET'
  402.      ('NOMC' (VNZ * RN) 'UZ' 'NATU' 'DISCRET' ) ;
  403.  
  404. RYN = RYH2 'ET' RYH2O 'ET' RYO2 ;
  405.  
  406. *
  407. * Calcul de l'energie thermique et du gamma
  408. *
  409.  
  410. RETHER = (A0H2 '*'  RN '*' TN) ;
  411.  
  412. LISCOM2 = 'MOTS' 'UX  ' 'UY  ' 'UZ  ';
  413. RECIN = 'PSCAL' (0.5 '*' VN) GN LISCOM2 LISCOM2 ;
  414.  
  415. RET = RECIN '+' RETHER ;
  416.  
  417.  CHUN  = 'MANUEL' 'CHPO' MAIL 1 'SCAL' 1.0 ;
  418.  
  419.  CV = 'NOMC' (A0H2 '*' CHUN) 'SCAL' 'NATU' 'DISCRET' ;
  420.  CP = CV '+' (PGAZ . 'H2  ' . 'R') ;
  421.  
  422.  GAMN = CP '/' CV ;
  423.  
  424.  RTOT = CP '-' CV ;
  425.  
  426.  PN = (RTOT '*' TN) '*' RN ;
  427.  
  428. *
  429. **** PRIM
  430. *
  431.  
  432. VITESSE PRES TEMPNEW GAMNEW
  433.        = 'PRIM' 'PERFTEMP' PGAZ RN GN RET ;
  434.  
  435. *
  436. **** TEST
  437. *
  438.  
  439. ERRV = ('MAXIMUM' (VN '-' VITESSE) 'ABS')
  440.         '/' ('MAXIMUM' VN) ;
  441.  
  442. ERRP  = 'MAXIMUM' ((PRES '-' PN)  '/'  PN) 'ABS';
  443.  
  444. ERRT  = 'MAXIMUM' ((TEMPNEW '-' TN) '/' TN) 'ABS';  
  445.  
  446. ERRG = 'MAXIMUM' ((GAMN '-' GAMNEW) '/' GAMN) 'ABS'; 
  447.  
  448.  
  449. 'SI' (('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG  ) '>'
  450.           ERRMAX));
  451.      'MESSAGE' ('CHAINE' 'Erreur maximum');
  452.      'LISTE' ('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG
  453.       ));
  454.      'ERREUR'  5;
  455. 'FINSI' ;
  456.  
  457. * Avec une temperature d'essai
  458.  
  459. VITESSE PRES TEMPNEW2 GAMNEW
  460.        = 'PRIM' 'PERFTEMP' PGAZ RN GN RET TEMPNEW;
  461.  
  462. ERRT  = 'MAXIMUM' ((TEMPNEW '-' TEMPNEW2) '/' TEMPNEW) 'ABS';
  463.  
  464. 'SI' (ERRT > ERRMAX);
  465.      'MESSAGE' ('CHAINE' 'Erreur maximum' ERRT) ;
  466.      'ERREUR'  5;
  467. 'FINSI' ;
  468.  
  469. **** On teste aussi la difference entre l'energie initiale et celle la.
  470.  
  471. RETHER_N = (A0H2 '*'  RN '*' TEMPNEW) ;
  472.  
  473. ERRET = 'MAXIMUM' ((RETHER '-' RETHER_N) '/' RETHER) 'ABS' ;
  474.  
  475.  
  476. 'SI' (ERRET > ERRMAX);
  477.      'MESSAGE' ('CHAINE' 'Erreur maximum' ERRET) ;
  478.      'ERREUR'  5;
  479. 'FINSI' ;
  480.  
  481. 'FIN' ;
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  

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