Télécharger la plaquette Cast3M


fr en

Télécharger primtest3_3D.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  primtest3_3D.dgibi
  2. ***********************************************************
  3. **** APPROCHE VF "Cell-Centred Formulation" pour la    ****
  4. **** solution des                                      ****
  5. **** Equations d'Euler pour un gaz parfait             ****
  6. **** OPERATEUR  PRIM                                   ****
  7. **** Cas: gaz multiespece "calorically perfect"        ****
  8. **** Cas 3D                                            ****
  9. ****                                                   ****
  10. **** A. BECCANTINI DRN/DMT/SEMT/TTMF    MARS  1998     ****
  11. ***********************************************************
  12.  
  13. 'OPTION'  'DIME' 3 ;
  14. 'OPTION'  'ELEM' CUB8 ;
  15. 'OPTION'  'ECHO'  0;
  16. 'OPTION'  'TRAC' 'X' ;
  17.  
  18. *
  19. *** GRAPHIQUES
  20. *
  21.  
  22. * GRAPH = VRAI ;
  23. GRAPH = FAUX ;
  24.  
  25.  
  26.  
  27. ************
  28. * MAILLAGE *
  29. ************
  30.  
  31. NX =  10 ;
  32. NY =  2  ;
  33. NZ =  2 ;
  34.  
  35. L  = 1.0D0  ;
  36. DX = L '/' NX '/' 2.0D0 ;
  37. H  = NY '*' DX  ;
  38. P  = NZ '*' DX ;
  39. xD =  0.5D0 '*' L ;
  40. xG = -1.0D0 '*' xD  ;
  41. yH =  0.5D0 '*' H   ;
  42. yB =  -1.0D0 '*' yH ;
  43. zV = 0.5D0 '*' P    ;
  44. zR = -1.0D0 '*' zV  ;
  45.  
  46. A1 = xG yB zR ;
  47. A2 = 0.0D0 yB zR ;
  48. A7 = xG yB  zV ;
  49. A8 = 0.0D0 yB  zV ;
  50.  
  51.  
  52. BAS1  = A1 'DROI' NX A2 ;
  53. BAS3  = A7 'DROI' NX A8 ;
  54. BAS5  = A1 'DROI' NZ A7 ;
  55. BAS6  = A2 'DROI' NZ A8 ;
  56.  
  57.  
  58. S11 = 'DALL' BAS3 BAS6 BAS1 BAS5 ;
  59.  
  60.  
  61. DOM1   = S11 'VOLU' NZ 'TRANS' (0.0 H 0.0) ;
  62. DOM2   = S11 'VOLU' NZ 'TRANS' (0.0 (-1.0*H) 0.0) ;
  63.  
  64.  
  65. DOMTOT = DOM1 ET DOM2;
  66. 'ELIMINATION' DOMTOT 1D-6;
  67.  
  68. $DOMTOT = 'DOMA' DOMTOT ;
  69. $DOM1   = 'DOMA' DOM1    'INCL' $DOMTOT ;
  70. $DOM2   = 'DOMA' DOM2    'INCL' $DOMTOT ;
  71.  
  72. 'SI' GRAPH;
  73.    'TRACER' ($DOMTOT . 'MAILLAGE' ) 'TITRE' 'Maillage';
  74. 'FINSI' ;
  75.  
  76. ************************
  77. **** MODELE DU GAZ  ****
  78. ************************
  79.  
  80. NESP = 4;
  81.  
  82. *
  83. *** GAS: H_2, O_2, H_2O, N_2
  84. *
  85. *   CP, CV en J/Kg/K @ T = 3000
  86. *
  87.  
  88. PGAS = 'TABLE' ;
  89.  
  90. PGAS . 'CP' = 'TABLE' ;
  91. PGAS . 'CP' . 'H2  '  = .18729066D+05 ;
  92. PGAS . 'CP' . 'O2  '  = .11886820D+04 ;
  93. PGAS . 'CP' . 'H2O '  = .31209047D+04 ;
  94. PGAS . 'CP' . 'N2  '  = .12993995D+04 ;
  95.  
  96. PGAS . 'CV' = 'TABLE' ;
  97. PGAS . 'CV' . 'H2  '  = .14571861D+05 ;
  98. PGAS . 'CV' . 'O2  '  = .92885670D+03 ;
  99. PGAS . 'CV' . 'H2O '  = .26589930D+04 ;
  100. PGAS . 'CV' . 'N2  '  =  .10024563D+04;
  101.  
  102. *
  103. **** Especes qui sont dans les equations d'Euler
  104. *
  105.  
  106. PGAS . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;
  107.  
  108. *
  109. **** Espece qui n'y est pas
  110. *
  111.  
  112.  
  113. PGAS . 'ESPNEULE' = 'N2  ';
  114.  
  115.  
  116.  
  117. *
  118. ***********************
  119. **** Les CHPOINTs  ****
  120. ***********************
  121.  
  122.  
  123. RN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' 1.9D0;
  124. RN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' 2.9D0;
  125. RN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
  126.  
  127. GN1   = 'KCHT' $DOM1   'VECT' 'CENTRE' (4.9D0 7.9D0 12.3);
  128. GN2   = 'KCHT' $DOM2   'VECT' 'CENTRE' (0.9D0 11.7D0 1.54);
  129. GN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN1 'ET' GN2);
  130.  
  131. PN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' 1.19D2;
  132. PN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' 2.92D3;
  133. PN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
  134.  
  135.  
  136. *
  137. **** Les fractions massiques YN 'ET' les masses RYN
  138. *
  139. *    N.B. YN n'as pas le meme supporte geometrique
  140. *         que RN, PN; en general la numerotation
  141. *         des noeuds est diferente, mais ce n'est
  142. *         pas grave: 'PRIM' change les masse des
  143. *         especes RYN dans la subroutine QUEPOI!!!
  144. *
  145.  
  146. UNCH =  'KCHT'  $DOM1 'SCAL' 'CENTRE' 1.0D0 ;
  147.  
  148. Y1      =  'KCHT'  $DOM1 'SCAL' 'CENTRE' 0.1   ;
  149. RY1     =  'NOMC' 'H2' (RN1 * Y1) 'NATU' 'DISCRET';
  150. YH2     =  'NOMC' 'H2' Y1 'NATU' 'DISCRET';
  151. Y2      =  'KCHT'  $DOM1 'SCAL' 'CENTRE' 0.2   ;
  152. RY2     =  'NOMC'  'O2' (RN1 * Y2) 'NATU' 'DISCRET';
  153. YO2     =  'NOMC' 'O2' Y2 'NATU' 'DISCRET';
  154. Y3      =  'KCHT'  $DOM1 'SCAL' 'CENTRE' 0.3   ;
  155. RY3     =  'NOMC' 'H2O' (RN1 * Y3) 'NATU' 'DISCRET';
  156. YH2O    =  'NOMC' 'H2O' Y3 'NATU' 'DISCRET';
  157. YTOT    =   Y1 '+' Y2 '+' Y3 ;
  158. Y4      =  'KCHT'  $DOM1 'SCAL' 'CENTRE' (UNCH '-' YTOT);
  159. YN2     =  'NOMC' 'N2' Y4 'NATU' 'DISCRET';
  160.  
  161. YND1 = YH2 'ET' YO2 'ET' YH2O 'ET' YN2;
  162. RYND1 = RY2 'ET' RY3 'ET' RY1;
  163.  
  164.  
  165. UNCH =  'KCHT'  $DOM2 'SCAL' 'CENTRE' 1.0D0 ;
  166. Y1      =  'KCHT'  $DOM2 'SCAL' 'CENTRE' 0.1   ;
  167. RY1     =  'NOMC' 'H2' (RN2 * Y1) 'NATU' 'DISCRET';
  168. YH2     =  'NOMC' 'H2' Y1 'NATU' 'DISCRET';
  169. Y2      =  'KCHT'  $DOM2 'SCAL' 'CENTRE' 0.2   ;
  170. RY2     =  'NOMC'  'O2' (RN2 * Y2) 'NATU' 'DISCRET';
  171. YO2     =  'NOMC' 'O2' Y2 'NATU' 'DISCRET';
  172. Y3      =  'KCHT'  $DOM2 'SCAL' 'CENTRE' 0.3   ;
  173. RY3     =  'NOMC' 'H2O' (RN2 * Y3) 'NATU' 'DISCRET';
  174. YH2O    =  'NOMC' 'H2O' Y3 'NATU' 'DISCRET';
  175. YTOT    =   Y1 '+' Y2 '+' Y3 ;
  176. Y4      =  'KCHT'  $DOM2 'SCAL' 'CENTRE' (UNCH '-' YTOT);
  177. YN2     =  'NOMC' 'N2' Y4 'NATU' 'DISCRET';
  178.  
  179. YND2 = YH2 'ET' YO2 'ET' YH2O 'ET' YN2;
  180. RYND2 = RY3 'ET' RY2 'ET' RY1;
  181.  
  182. YN = YND1 'ET' YND2;
  183. RYN = RYND1 'ET' RYND2 ;
  184.  
  185. *
  186. **** Control de la Table PGAS
  187. *
  188.  
  189.  
  190. *
  191. **** La vitesse
  192. *
  193.  
  194. VITX = ('EXCO' 'UX  ' GN) '/' RN;
  195. VITY = ('EXCO' 'UY  ' GN) '/' RN;
  196. VITZ = ('EXCO' 'UZ  ' GN) '/' RN;
  197.  
  198. *
  199. **** Les GAMMAs
  200. *
  201.  
  202.  NOMCOM = 'EXTRAIRE' YN 'COMP';
  203.  CPTOT = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 0.0D0;
  204.  CVTOT = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 0.0D0;
  205. 'REPETER' BL1 ('DIME' NOMCOM);
  206.     NOMCEL = 'EXTRAIRE' NOMCOM &BL1;
  207.     YCEL   = 'EXCO' NOMCEL YN ;
  208.     CPCEL  = (PGAS . 'CP' . NOMCEL);
  209.     CVCEL  = (PGAS . 'CV' . NOMCEL);
  210.     CPTOT = CPTOT '+' (YCEL * CPCEL);
  211.     CVTOT = CVTOT '+' (YCEL * CVCEL);
  212. 'FIN' BL1;
  213.  GAMMA = CPTOT '/' CVTOT;
  214.  RGAS = (CPTOT '-' CVTOT);
  215.  
  216.  
  217. *
  218. *** L'energie totale (ROET) et  la temperature
  219. *
  220.  
  221. GM1 = GAMMA '-' ('KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.0D0);
  222. ETHER = PN '/' GM1;
  223. lcg=mots 'UX' 'UY' 'UZ';
  224. ECIN = (0.5D0 * (GN lcg 'PSCA' GN lcg)) '/' RN;
  225. EN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (ECIN '+' ETHER);
  226. TEMPE = (PN '/' RN) '/' RGAS;
  227.  
  228.  
  229. ************************
  230. ****  L'operateur  *****
  231. ************************
  232.  
  233. VITESSE PRES TEMPNEW YNNEW GAMNEW
  234.        = 'PRIM' 'PERFMULT' PGAS RN GN EN RYN ;
  235.  
  236.  
  237. *
  238. **** L'erreur (pas d'integral pour semplifier le calcul);
  239. *
  240.  
  241. VIT1X = 'EXCO' 'UX  ' VITESSE ;
  242. VIT1Y = 'EXCO' 'UY  ' VITESSE ;
  243. VIT1Z = 'EXCO' 'UZ  ' VITESSE ;
  244.  
  245. ERRVX = 'MAXIMUM' (VIT1X '-' VITX) 'ABS' ;
  246. ERRVY = 'MAXIMUM' (VIT1Y '-' VITY) 'ABS' ;
  247. ERRVZ = 'MAXIMUM' (VIT1Z '-' VITZ) 'ABS';
  248. VCELL = ('MAXIMUM' VITX 'ABS') '+'
  249.         ('MAXIMUM' VITY 'ABS' )'+'
  250.         ('MAXIMUM' VITZ 'ABS' );
  251. 'SI' (VCELL > 0);
  252.    ERRVX = ERRVX '/' VCELL;
  253.    ERRVY = ERRVY '/' VCELL;
  254.    ERRVZ = ERRVZ '/' VCELL;
  255. 'FINSI'  ;
  256.  
  257. ERRP  = ('MAXIMUM' (PRES '-' PN) 'ABS') '/' ('MAXIMUM' PN);
  258.  
  259. ERRT  = ('MAXIMUM' (TEMPNEW '-' TEMPE) 'ABS') '/'
  260.          ('MAXIMUM' TEMPE) ;
  261.  
  262. ERRG = ('MAXIMUM' (GAMMA '-' GAMNEW) 'ABS') '/'
  263.        ('MAXIMUM' GAMMA);
  264.  
  265. LCOMP =  'EXTRAIRE' YNNEW 'COMP' ;
  266. NCOMP = 'DIME' LCOMP ;
  267.  
  268. ERRY = 0.0D0 ;
  269. 'REPETER' BL1 NCOMP ;
  270.       NOMCEL = 'EXTRAIRE' LCOMP &BL1 ;
  271.       YN1 = 'EXCO' NOMCEL YN 'SCAL'  ;
  272.       YN2 = 'EXCO' NOMCEL YNNEW 'SCAL' ;
  273.       ERRY  = ERRY '+' ('MAXIMUM' (YN1 '-' YN2) 'ABS') ;
  274. 'FINSI' ;
  275.  
  276. ERRY = ERRY '/' NCOMP;
  277.  
  278.  
  279.  
  280. 'SI' (('MAXIMUM' ('PROG' ERRVX ERRVY ERRVZ ERRP ERRT ERRG ERRY ) '>'
  281.           1.0D-12));
  282.      'MESSAGE' ('CHAINE' 'Erreur maximum');
  283.      'LISTE' ('MAXIMUM' ('PROG' ERRVX ERRVY ERRVZ ERRP ERRT ERRG ERRY
  284.      ));
  285.     'ERREUR'  5;
  286. 'FINSI' ;
  287.  
  288. 'FIN' ;
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  

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