Télécharger pret_ther2.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : pret_ther2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. **********************************************************************
  5. **** APPROCHE VF "Cell-Centred Formulation" pour la solution des ****
  6. **** Equations d'Euler pour un gaz parfait. ****
  7. **** ****
  8. **** OPERATEUR PRET ****
  9. **** Operateur qui 'recontruit les variables primitives aux faces ****
  10. **** Cas gaz "thermally perfect" multi-especes ****
  11. **** Premier ordre en espace, premier ordre en temps ****
  12. **** Enterieur et murs ****
  13. **** ****
  14. **** A. BECCANTINI DRN/DMT/SEMT/TTMF FEVRIER 2000 ****
  15. **********************************************************************
  16.  
  17. 'OPTION' 'DIME' 2 ;
  18. 'OPTION' 'ELEM' QUA4 ;
  19. 'OPTION' 'ECHO' 0 ;
  20. 'OPTION' 'TRAC' 'PS' ;
  21.  
  22. *
  23. *** GRAPH
  24. *
  25.  
  26. * GRAPH = VRAI ;
  27. GRAPH = FAUX ;
  28.  
  29. ***************************
  30. *** PROCEDURE RUO1 ****
  31. ***************************
  32. *
  33. *** Precedure pour la rotation d'un tenseur de premier
  34. * ordre
  35. *
  36.  
  37. 'DEBPROC' RUO1 ;
  38. 'ARGUMENT' ALPHA*'FLOTTANT' UN*'FLOTTANT' UT*'FLOTTANT';
  39.  
  40. SINA = 'SIN' ALPHA ;
  41. COSA = 'COS' ALPHA ;
  42.  
  43. UX = (UN * COSA ) '-' (UT * SINA) ;
  44. UY = (UN * SINA ) '+' (UT * COSA) ;
  45.  
  46.  
  47. 'FINPROC' UX UY;
  48.  
  49. ***************************
  50. *** PROCEDURE RUO2 ****
  51. ***************************
  52. *
  53. *** Precedure pour la rotation d'un tenseur de deuxieme
  54. * ordre
  55. *
  56.  
  57. 'DEBPROC' RUO2 ;
  58. 'ARGUMENT' ALPHA*'FLOTTANT' UNN*'FLOTTANT' UNT*'FLOTTANT'
  59. UTN*'FLOTTANT' UTT*'FLOTTANT';
  60. *
  61. **** (n,t) -> (x,y)
  62. *
  63. * n = CA x '+' SA y ;
  64. * t = -SA x '+' CA y ;
  65. *
  66. * UNT = DUN/DT
  67. *
  68.  
  69. SA = 'SIN' ALPHA ;
  70. CA = 'COS' ALPHA ;
  71. CA2 = CA * CA ;
  72. CASA = CA * SA ;
  73. SA2 = SA * SA;
  74.  
  75. UXX = (CA2 * UNN) '-' (CASA * (UNT '+' UTN)) '+' (SA2 *UTT) ;
  76. UYX = (CASA * (UNN '-' UTT)) '-' (SA2 * UNT) '+' (CA2 * UTN ) ;
  77. UXY = (CASA * (UNN '-' UTT)) '+' (CA2 * UNT) '-' (SA2 * UTN ) ;
  78. UYY = (SA2 * UNN) '+' (CASA * (UNT '+' UTN)) '+' (CA2 *UTT) ;
  79.  
  80. 'FINPROC' UXX UXY UYX UYY;
  81.  
  82. ***************************
  83. ***** DOMAINE SPATIAL ****
  84. ***************************
  85.  
  86. *
  87. **** Deux carre
  88. *
  89.  
  90. A1 = 0.0D0 0.0D0;
  91. A2 = 1.0D0 0.0D0;
  92. A3 = 2.0D0 0.0D0;
  93. A4 = 2.0D0 1.0D0;
  94. A5 = 1.0D0 1.0D0;
  95. A6 = 0.0D0 1.0D0;
  96.  
  97. L12 = A1 'DROIT' 1 A2;
  98. L23 = A2 'DROIT' 1 A3;
  99. L34 = A3 'DROIT' 1 A4;
  100. L45 = A4 'DROIT' 1 A5;
  101. L56 = A5 'DROIT' 1 A6;
  102. L61 = A6 'DROIT' 1 A1;
  103. L25 = A2 'DROIT' 1 A5;
  104.  
  105.  
  106. DOM10 = 'DALL' L12 L25 L56 L61
  107. 'PLANE';
  108. DOM20 = 'DALL' L23 L34 L45 ('INVERSE' L25)
  109. 'PLANE';
  110.  
  111. *
  112. *** Point face entre le deux carre, ou on fait les controlles
  113. *
  114.  
  115. P10 = 1.0 0.5 ;
  116.  
  117. *
  118. *** Point face sur le mur (a droite)
  119. *
  120.  
  121. P20 = 2.0 0.5 ;
  122.  
  123.  
  124. *
  125. *** Etats a gauche et a droite du point face P10
  126. *
  127. * ro, un, ut, p
  128.  
  129.  
  130.  
  131. rog = 1.100 ;
  132. pg = 2.13D0 ;
  133. ung = 1.17D0 ;
  134. utg = 1.14D0 ;
  135. yh2g = 0.21 ;
  136. sca1g = 1.76 ;
  137. sca2g = 311.76 ;
  138.  
  139.  
  140.  
  141. rod = 1.50 ;
  142. pd = 2.5D0 ;
  143. und = 3.5D0 ;
  144. utd = 4.5D0 ;
  145. yh2d = 0.23 ;
  146. sca1d = 217.6 ;
  147. sca2d = 110.0 ;
  148.  
  149. ETATG = 'PROG' rog ung utg pg yh2g sca1g sca2g ;
  150. ETATD = 'PROG' rod und utd pd yh2d sca1d sca2d ;
  151.  
  152. **** Cas multiespece : la table proprieté de gaz
  153. *
  154. * PGAZ2 avec 'SCALPASS'
  155.  
  156. PGAZ = 'TABLE' ;
  157. PGAZ2 = 'TABLE' ;
  158.  
  159. *
  160. **** Especes qui sont dans les equations d'Euler
  161. *
  162.  
  163. PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' ;
  164. PGAZ2 . 'ESPEULE' = 'MOTS' 'H2 ' ;
  165.  
  166. *
  167. **** Espece qui n'y est pas
  168. *
  169.  
  170.  
  171. PGAZ . 'ESPNEULE' = 'O2 ';
  172. PGAZ2 . 'ESPNEULE' = 'O2 ';
  173.  
  174.  
  175. *
  176.  
  177. PGAZ . 'H2 ' = 'TABLE' ;
  178. PGAZ . 'O2 ' = 'TABLE' ;
  179. PGAZ2 . 'H2 ' = 'TABLE' ;
  180. PGAZ2 . 'O2 ' = 'TABLE' ;
  181.  
  182. *
  183. **** R (J/Kg/K)
  184. *
  185.  
  186. PGAZ . 'H2 ' . 'R' = 4130.0 ;
  187. PGAZ . 'O2 ' . 'R' = 259.8 ;
  188. PGAZ2 . 'H2 ' . 'R' = 4130.0 ;
  189. PGAZ2 . 'O2 ' . 'R' = 259.8 ;
  190.  
  191. *
  192. **** Regressions polynomials
  193. *
  194.  
  195. PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.0 0.0 0.0 0.0 ;
  196. PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.0 0.0 0.0 0. ;
  197. PGAZ2 . 'H2 ' . 'A' = 'PROG' 9834.91866 0.0 0.0 0.0 0.0 ;
  198. PGAZ2 . 'O2 ' . 'A' = 'PROG' 575.012333 0.0 0.0 0.0 0. ;
  199.  
  200. *
  201. **** "Enthalpies" (ou energies) de formations a OK (J/Kg)
  202. * Note: ce sont des entites numeriques
  203. * h_i = h_i(T0) '-' \int_0^{T0} cp_i(x) dx =
  204. * h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
  205. *
  206.  
  207. PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
  208. PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
  209. PGAZ2 . 'H2 ' . 'H0K' = -4.195D6 ;
  210. PGAZ2 . 'O2 ' . 'H0K' = -2.634D5 ;
  211.  
  212. PGAZ2 . 'SCALPASS' = 'MOTS' 'TITU' 'TOTO' ;
  213.  
  214.  
  215. **** Test1: rotation
  216.  
  217. ANGLE = 7.0D0;
  218. DANGLE = 85;
  219. ORIG = 0.0D0 0.0D0;
  220.  
  221.  
  222. 'REPETER' BL1 5;
  223. ANGLE = ANGLE '+' DANGLE ;
  224.  
  225. uxg uyg = RUO1 ANGLE UNG UTG ;
  226. uxd uyd = RUO1 ANGLE UND UTD ;
  227.  
  228. 'MESSAGE' ;
  229. 'MESSAGE' (CHAIN 'Angle de rotation= ' ANGLE);
  230. 'MESSAGE' ;
  231.  
  232. DOM1 = DOM10 'TOURNER' ANGLE ORIG;
  233. DOM2 = DOM20 'TOURNER' ANGLE ORIG;
  234. P1 = P10 'TOURNER' ANGLE ORIG;
  235. P2 = P20 'TOURNER' ANGLE ORIG;
  236.  
  237. DOMTOT = DOM1 ET DOM2;
  238. 'ELIMINATION' (DOM1 ET DOM2) 1D-6;
  239. DOMTOT = DOM1 ET DOM2;
  240.  
  241.  
  242. $DOMTOT = 'MODELISER' DOMTOT 'EULER';
  243.  
  244. $DOM1 = 'MODELISER' DOM1 'EULER';
  245. $DOM2 = 'MODELISER' DOM2 'EULER';
  246.  
  247. TDOMTOT = 'DOMA' $DOMTOT 'VF';
  248.  
  249. TDOM1 = 'DOMA' $DOM1 'VF';
  250. TDOM2 = 'DOMA' $DOM2 'VF';
  251.  
  252. MDOM1 = TDOM1 . 'QUAF' ;
  253. MDOM2 = TDOM2 . 'QUAF' ;
  254. MDOMTOT = TDOMTOT . 'QUAF' ;
  255.  
  256. 'ELIMINATION' (MDOMTOT ET MDOM1) 0.0001 ;
  257. 'ELIMINATION' (MDOMTOT ET MDOM2) 0.0001 ;
  258.  
  259. 'SI' GRAPH;
  260. 'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL')
  261. 'ET' P10) 'TITRE' 'Domaine et FACEL';
  262. 'FINSI' ;
  263.  
  264. **** La densité, la pression, la vitesse
  265.  
  266. RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog ;
  267. RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod ;
  268. RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2) ;
  269.  
  270. PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pg ;
  271. PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pd ;
  272. PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2) ;
  273.  
  274. VN1X = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'UX ' uxg ;
  275. VN1Y = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'UY ' uyg ;
  276. VN2X = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'UX ' uxd ;
  277. VN2Y = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'UY ' uyd ;
  278. VN = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
  279. (VN1X 'ET' VN1Y 'ET' VN2X 'ET' VN2Y);
  280.  
  281. YN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2 ' yh2g ;
  282. YN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'H2 ' yh2d ;
  283. YN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2 ' (YN1 'ET' YN2) ;
  284.  
  285. SCAN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP'
  286. ('MOTS' 'TITU' 'TOTO') (sca1g sca2g) ;
  287. SCAN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP'
  288. ('MOTS' 'TITU' 'TOTO') (sca1d sca2d) ;
  289. SCAN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP'
  290. ('MOTS' 'TITU' 'TOTO') (SCAN1 'ET' SCAN2) ;
  291.  
  292. *
  293. **** Cas monoespece: dans la table PGAZ il n'y a pas de
  294. * PGAZ . 'ESPEULE'
  295.  
  296.  
  297. ROF VITF PF YF = 'PRET' 'PERFTEMP' 1 1 $DOMTOT PGAZ RN VN PN YN ;
  298.  
  299.  
  300. ROF2 VITF2 PF2 YF2 SCAF = 'PRET' 'PERFTEMP'
  301. 1 1 $DOMTOT PGAZ2 RN VN PN YN SCAN ;
  302.  
  303. ERRO = 'MAXIMUM' (ROF '-' ROF2) 'ABS' ;
  304. ERVI = 'MAXIMUM' (VITF '-' VITF2) 'ABS' ;
  305. ERP = 'MAXIMUM' (PF '-' PF2) 'ABS' ;
  306. ERY = 'MAXIMUM' (YF '-' YF2) 'ABS' ;
  307.  
  308.  
  309. 'SI' (('MAXIMUM' ('PROG' ERRO ERVI ERP ERY)) > 1.0D-10) ;
  310. 'MESSAGE' 'Erreur trop importante' ;
  311. 'ERREUR' 5 ;
  312. 'FINSI' ;
  313.  
  314. *
  315. **** Redefinition de P1, P2 dans $DOMTOT 'FACE'
  316. *
  317.  
  318. P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
  319. P2 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P2;
  320.  
  321. *********************************************************
  322. *** Control des etats sur la surface qui contient P1 ****
  323. *********************************************************
  324.  
  325.  
  326. GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
  327. GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
  328.  
  329. ROGEOP1 = 'REDU' ROF GEOP1 ;
  330. VGEOP1 = 'REDU' VITF GEOP1 ;
  331. PGEOP1 = 'REDU' PF GEOP1 ;
  332. REFGEOP1 = 'REDU' VITF GEOP2 ;
  333. YGEOP1 = 'REDU' YF GEOP1 ;
  334. SGEOP1 = 'REDU' SCAF GEOP1 ;
  335.  
  336. rofag = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 1;
  337. rofad = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 3;
  338. unfag = 'EXTRAIRE' VGEOP1 'UN ' 1 1 1;
  339. unfad = 'EXTRAIRE' VGEOP1 'UN ' 1 1 3;
  340. utfag = 'EXTRAIRE' VGEOP1 'UT ' 1 1 1;
  341. utfad = 'EXTRAIRE' VGEOP1 'UT ' 1 1 3;
  342. pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
  343. pfad = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 3;
  344. yfag = 'EXTRAIRE' YGEOP1 'H2 ' 1 1 1;
  345. yfad = 'EXTRAIRE' YGEOP1 'H2 ' 1 1 3;
  346. s1fag = 'EXTRAIRE' SGEOP1 'TITU' 1 1 1;
  347. s1fad = 'EXTRAIRE' SGEOP1 'TITU' 1 1 3;
  348. s2fag = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 1;
  349. s2fad = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 3;
  350.  
  351. *
  352. *
  353. **** Orientation de la normal n de castem par raport a la
  354. * notre; t est par consequence
  355. *
  356.  
  357. NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1;
  358. NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1;
  359.  
  360. 'SI' (('ABS' NCOS) > ('ABS' NSIN));
  361. ORIENT = ('COS' ANGLE) '/' NCOS;
  362. 'SINON';
  363. ORIENT = ('SIN' ANGLE) '/' NSIN;
  364. 'FINSI' ;
  365.  
  366. ORIENT = 'SIGNE' 'FLOTTANT' ORIENT;
  367.  
  368. *
  369. **** ORIENT = -1 -> Mon etat gauche est son etat droite
  370. *
  371.  
  372. 'SI' (ORIENT > 0);
  373.  
  374. ERRLIG = 'PROG' rofag (unfag '*' ORIENT)
  375. (utfag '*' ORIENT) pfag yfag s1fag s2fag ;
  376.  
  377. ERRLID = 'PROG' rofad (unfad '*' ORIENT)
  378. (utfad '*' ORIENT) pfad yfad s1fad s2fad ;
  379.  
  380. 'SINON' ;
  381.  
  382. ERRLID = 'PROG' rofag (unfag '*' ORIENT)
  383. (utfag '*' ORIENT) pfag yfag s1fag s2fag ;
  384.  
  385. ERRLIG = 'PROG' rofad (unfad '*' ORIENT)
  386. (utfad '*' ORIENT) pfad yfad s1fad s2fad ;
  387.  
  388. 'FINSI' ;
  389.  
  390. ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ETATG '-' ERRLIG) 'ABS')
  391. ('MAXIMUM' (ETATD '-' ERRLID) 'ABS'));
  392.  
  393.  
  394. 'SI' (ERRO > 1.0D-12) ;
  395. 'MESSAGE' 'Ordre en espace = 1';
  396. 'MESSAGE' 'Ordre en temps = 1';
  397. 'ERREUR' 5 ;
  398. 'FINSI' ;
  399.  
  400. *********************************************************
  401. *** Control des etats sur la surface qui contient P2 ****
  402. *********************************************************
  403.  
  404.  
  405. GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
  406. GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
  407.  
  408.  
  409.  
  410. ROGEOP1 = 'REDU' ROF GEOP1;
  411. VGEOP1 = 'REDU' VITF GEOP1;
  412. PGEOP1 = 'REDU' PF GEOP1;
  413. REFGEOP1 = 'REDU' VITF GEOP2;
  414. YGEOP1 = 'REDU' YF GEOP1 ;
  415. SGEOP1 = 'REDU' SCAF GEOP1 ;
  416.  
  417. rofag = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 1;
  418. rofad = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 3;
  419. unfag = 'EXTRAIRE' VGEOP1 'UN ' 1 1 1;
  420. unfad = 'EXTRAIRE' VGEOP1 'UN ' 1 1 3;
  421. utfag = 'EXTRAIRE' VGEOP1 'UT ' 1 1 1;
  422. utfad = 'EXTRAIRE' VGEOP1 'UT ' 1 1 3;
  423. pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
  424. pfad = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 3;
  425. yfag = 'EXTRAIRE' YGEOP1 'H2 ' 1 1 1;
  426. yfad = 'EXTRAIRE' YGEOP1 'H2 ' 1 1 3;
  427. s1fag = 'EXTRAIRE' SGEOP1 'TITU' 1 1 1;
  428. s1fad = 'EXTRAIRE' SGEOP1 'TITU' 1 1 3;
  429. s2fag = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 1;
  430. s2fad = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 3;
  431.  
  432. * Par convention l'etat droite est l'etat mirroir
  433. *
  434. * Sur le mur: Gauche = etat reel
  435. * Droit = etat mirroir
  436. *
  437.  
  438. LOG1 = ('ABS' (rod '-' rofag)) > 1.0D-12 ;
  439. LOG1 = 'OU' LOG1 (('ABS' (rod '-' rofad)) > 1.0D-12) ;
  440. LOG1 = 'OU' LOG1 (('ABS' (und '-' unfag)) > 1.0D-12) ;
  441. LOG1 = 'OU' LOG1 (('ABS' (und '+' unfad)) > 1.0D-12) ;
  442. LOG1 = 'OU' LOG1 (('ABS' (utd '-' utfag)) > 1.0D-12) ;
  443. LOG1 = 'OU' LOG1 (('ABS' (utd '-' utfad)) > 1.0D-12) ;
  444. LOG1 = 'OU' LOG1 (('ABS' (pd '-' pfag)) > 1.0D-12) ;
  445. LOG1 = 'OU' LOG1 (('ABS' (pd '-' pfad)) > 1.0D-12) ;
  446. LOG1 = 'OU' LOG1 (('ABS' (yh2d '-' yfag)) > 1.0D-12) ;
  447. LOG1 = 'OU' LOG1 (('ABS' (yh2d '-' yfad)) > 1.0D-12) ;
  448. LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fag)) > 1.0D-12) ;
  449. LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fad)) > 1.0D-12) ;
  450. LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fag)) > 1.0D-12) ;
  451. LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fad)) > 1.0D-12) ;
  452. LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fag)) > 1.0D-12) ;
  453. LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fad)) > 1.0D-12) ;
  454. LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fag)) > 1.0D-12) ;
  455. LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fad)) > 1.0D-12) ;
  456.  
  457.  
  458. 'SI' (LOG1) ;
  459. 'MESSAGE' 'Mur' ;
  460. 'MESSAGE' 'Ordre en espace = 1';
  461. 'MESSAGE' 'Ordre en temps = 1';
  462. 'ERREUR' 5 ;
  463. 'FINSI' ;
  464.  
  465. 'FIN' BL1 ;
  466.  
  467.  
  468.  
  469. 'FIN' ;
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  

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