Télécharger pret_ther.dgibi

Retour à la liste

Numérotation des lignes :

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

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