Télécharger la plaquette Cast3M


fr en

Télécharger konv_ther_cons.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  konv_ther_cons.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***********************************************************
  5. **** APPROCHE VF "Cell-Centred Formulation" pour la    ****
  6. **** solution des                                      **** 
  7. **** Equations d'Euler pour un gaz parfait             ****
  8. **** OPERATEURS PRET, KONV                             ****
  9. ****                                                   ****
  10. **** Cas gaz monoespece  "thermally perfect"           ****
  11. **** Transport de scalaires passifs                    ****
  12. **** Consistence                                       ****
  13. ****                                                   ****
  14. **** Methodes: VLH, Colella-Glaz                       ****
  15. ****                                                   ****
  16. **** A. BECCANTINI DRN/DMT/SEMT/TTMF    FEVRIER 2000   ****
  17. ***********************************************************
  18.  
  19. 'OPTION'  'DIME' 2 ;
  20. 'OPTION'  'ELEM' QUA4 ;
  21. 'OPTION'  'ECHO'  0 ;
  22. 'OPTION'  'TRAC' 'X' ;
  23.  
  24. *
  25. *** GRAPH
  26. *
  27.  
  28. GRAPH = FAUX ;
  29. * GRAPH = VRAI ;
  30.  
  31. **** Cas monoespece : la table proprieté de gaz
  32.  
  33. PGAZ = 'TABLE' ;
  34. PGAZ2 = 'TABLE' ;
  35.  
  36. *
  37. **** Ordre des polynoms cv_i
  38. *
  39.  
  40. PGAZ . 'NORD' = 4 ;
  41. PGAZ2 . 'NORD' = 4 ;
  42.  
  43. *
  44. **** Especes qui sont dans les equations d'Euler (rien)
  45. *
  46.  
  47. *
  48. **** Espece qui n'y est pas
  49. *
  50.  
  51.  
  52. PGAZ . 'ESPNEULE' = 'H2  ';
  53. PGAZ2 . 'ESPNEULE' = 'H2  ';
  54.  
  55. *
  56.  
  57. PGAZ .  'H2  ' = 'TABLE'  ;
  58. PGAZ2 .  'H2  ' = 'TABLE'  ;
  59.  
  60. *
  61. **** R (J/Kg/K)
  62. *
  63.  
  64. PGAZ .  'H2  ' . 'R' = 4130.0 ;
  65. PGAZ2 .  'H2  ' . 'R' = 4130.0 ;
  66.  
  67. *
  68. **** Regressions polynomials
  69. *
  70.  
  71. PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
  72.                                -2.37281455E-07 1.84701105E-11 ;
  73. PGAZ2 .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
  74.                                -2.37281455E-07 1.84701105E-11 ;
  75.  
  76. *
  77. **** "Enthalpies" (ou energies) de formations a OK (J/Kg)
  78. *     Note: ce sont des entites numeriques
  79. *
  80.  
  81. PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
  82. PGAZ2 .  'H2  ' . 'H0K' = -4.195D6 ;
  83.  
  84.  
  85. * PGAZ2: transport de scalaires passifs
  86.  
  87. PGAZ2 . 'SCALPASS' = 'MOTS' 'TUTU' 'TATA' ;
  88.  
  89. *************************** 
  90. ***** DOMAINE SPATIAL  ****
  91. ***************************
  92.  
  93.  
  94. A1 = 0.0D0 0.0D0;
  95. A2 = 1.0D0 0.0D0;
  96. A3 = 2.0D0 0.0D0;
  97. A4 = 2.0D0 1.0D0;
  98. A5 = 1.0D0 1.0D0;
  99. A6 = 0.0D0 1.0D0;
  100.  
  101. L12 = A1 'DROIT' 1 A2;
  102. L23 = A2 'DROIT' 1 A3;
  103. L34 = A3 'DROIT' 1 A4;
  104. L45 = A4 'DROIT' 1 A5;
  105. L56 = A5 'DROIT' 1 A6;
  106. L61 = A6 'DROIT' 1 A1;
  107. L25 = A2 'DROIT' 1 A5;
  108.  
  109.  
  110. DOM10  = 'DALL' L12 L25 L56 L61
  111.         'PLANE';
  112. DOM20  = 'DALL' L23 L34 L45 ('INVERSE' L25)
  113.         'PLANE';
  114. *
  115. *** Point ou on controlle la consistence
  116. *
  117.  
  118. P10 = 1.0 0.5;
  119.  
  120. *
  121. *** Etats gauche et droite 
  122. *
  123.  
  124. rogd = 1.412 ;
  125. pgd = 1.17D6 ;
  126. ungd  = 8.34 ;
  127. utgd  = 18.2 ;
  128. tutugd = 14.9 ;
  129. tatagd = 23.8 ;
  130.  
  131. *
  132. **** tg
  133. *
  134.  
  135. tgd = pgd '/' (rogd * (PGAZ .  'H2  ' . 'R')) ;
  136.  
  137.  
  138. *
  139. **** gamg, gamd , htg , htd
  140. *
  141. *
  142.  
  143. A0H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 1 ;
  144. A1H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 2 ;
  145. A2H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 3 ;
  146. A3H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 4 ;
  147. A4H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 5 ;
  148.  
  149. T = tgd ;
  150. T2 = T * T ;
  151. T3 = T2 * T ;
  152. T4 = T3 * T ;
  153. T5 = T4 * T ;
  154.  
  155.  cvgd =  A0H2 '+' (A1H2 * T) '+' (A2H2 * T2) '+' (A3H2 * T3) '+'
  156.         (A4H2 * T4) ;
  157.  
  158.  ETGd = (A0H2 * T) '+' (A1H2 * T2 '/' 2) '+' (A2H2 * T3 '/' 3) '+'
  159.       (A3H2 * T4 '/' 4) '+' (A4H2 * T5 '/' 5) ;
  160.  
  161.  ecingd = 0.5D0 '*'  ((ungd '*' ungd) '+' (utgd '*' utgd));
  162.  
  163.  retgd = rogd '*' (ETGd '+' ecingd) ;
  164.  
  165.  gamgd = (cvgd '+' (PGAZ .  'H2  ' . 'R')) '/' cvgd ;
  166.  
  167. *
  168. ***  flux en (n,t)
  169. *
  170.  
  171. f1gd = ungd  '*'  rogd ;
  172. f2gd = (rogd '*' (ungd '*' ungd)) '+' pgd ;
  173. f3gd = rogd '*' ungd '*' utgd ;
  174. f4gd = ungd '*' (retgd '+' pgd);
  175. ftutugd =  ungd  '*'  rogd '*' tutugd ;
  176. ftatagd =  ungd  '*'  rogd '*' tatagd ;
  177.  
  178.  
  179.  
  180. ****************************************************
  181. ****************************************************
  182. ******** Boucle sur les angles             *********
  183. ****************************************************
  184. ****************************************************
  185.  
  186. DANGLE = 360 '/' 7.15;
  187. ANGLE = 11.3 ;
  188.  
  189. 'REPETER' BLOC 8;
  190.  
  191. *
  192. *** Rotation
  193. *
  194.  
  195.  
  196. ANGLE = ANGLE '+' DANGLE;
  197. ORIG = 0.0D0 0.0D0;
  198.  
  199. MESSAGE;
  200. MESSAGE (CHAIN 'Angle de rotation= ' ANGLE);
  201. MESSAGE;
  202.  
  203. DOM1   = DOM10   'TOURNER' ANGLE ORIG;
  204. DOM2   = DOM20   'TOURNER' ANGLE ORIG;
  205. P1     = P10     'TOURNER' ANGLE ORIG;
  206.  
  207.     'ELIMINATION' (DOM1 ET DOM2) 1D-6;
  208.     DOMTOT = DOM1 ET DOM2;
  209.  
  210.  $DOMTOT = 'MODELISER' DOMTOT 'EULER';
  211.  $DOM1 = 'MODELISER'  DOM1 'EULER';
  212.  $DOM2 = 'MODELISER'  DOM2 'EULER';
  213.  TDOMTOT = 'DOMA'  $DOMTOT 'VF';
  214.  TDOM1 = 'DOMA'  $DOM1 'VF';
  215.  TDOM2 = 'DOMA'  $DOM2 'VF';
  216.  MDOM1 = TDOM1 . 'QUAF' ;
  217.  MDOM2 = TDOM2 . 'QUAF' ;
  218.  MDOMTOT = TDOMTOT . 'QUAF' ;
  219.  
  220.  'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) 1.E-6 ;
  221.  
  222.  'SI' GRAPH;
  223.       'TRACER' (('DOMA' $DOMTOT 'MAILLAGE' ) 'ET' 
  224.       ('DOMA' $DOMTOT 'FACEL') 'ET' P1) 'TITRE' 'Domaine et FACEL';
  225.  'FINSI' ;
  226.  
  227. *
  228. **** Redefinition de P1 dans $DOMTOT 'FACE'
  229. *
  230.  
  231. P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
  232.  
  233. ***********************
  234. **** Les CHPOINTs  ****
  235. ***********************
  236.  
  237. uxgd = (ungd '*' ('COS' ANGLE)) '-' (utgd '*' ('SIN' ANGLE));
  238. uygd = (ungd '*' ('SIN' ANGLE)) '+' (utgd '*' ('COS' ANGLE));
  239.  
  240.  
  241. RN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' rogd;
  242. RN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' rogd;;
  243. RN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
  244.  
  245. VN1   = 'KCHT' $DOM1   'VECT' 'CENTRE' (uxgd uygd);
  246. VN2   = 'KCHT' $DOM2   'VECT' 'CENTRE' (uxgd uygd);
  247. VN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VN1 'ET' VN2);
  248.  
  249. PN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' pgd;
  250. PN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' pgd;
  251. PN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
  252.  
  253. SCAN1   = 'KCHT' $DOM1   'VECT' 'CENTRE'
  254.           'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ;
  255. SCAN2   = 'KCHT' $DOM2   'VECT' 'CENTRE'
  256.           'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ;
  257. SCAN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
  258.           'COMP' (PGAZ2 . 'SCALPASS') (SCAN1 'ET' SCAN2);
  259.  
  260. ***************************
  261. ****  L'operateur  PRET****
  262. ***************************
  263.  
  264. ORDESP = 1;
  265. ORDTEM = 1;
  266.  
  267.  
  268. ROF VITF PF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
  269.     $DOMTOT PGAZ RN VN PN ;
  270.  
  271. ROF2 VITF2 PF2 SCAF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
  272.     $DOMTOT PGAZ2 RN VN PN SCAN ;
  273.  
  274. GEOP1    = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
  275. GEOP2    = ('DOMA' $DOMTOT 'FACE')  'ELEM' 'APPUYE' 'LARGEMENT' P1;
  276.  
  277.  
  278. REFGEOP1 = 'REDU' VITF GEOP2;
  279.  
  280. *
  281. **** Orientation de la normal n de castem par raport a la
  282. *    notre; t est par consequence
  283. *
  284.  
  285. NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1;
  286. NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1;
  287.  
  288. 'SI' (('ABS' NCOS) > ('ABS' NSIN));
  289.    ORIENT = ('COS' ANGLE) '/' NCOS;
  290. 'SINON'; 
  291.    ORIENT = ('SIN' ANGLE) '/' NSIN;
  292. 'FINSI' ;
  293.  
  294. ORIENT = 'SIGN' ORIENT;   
  295.  
  296. **************************
  297. **** 'KONV', TEST VLH ****
  298. **************************
  299.  
  300. TFLUX DT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' 'VLH'
  301.   $DOMTOT PGAZ ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')  ROF VITF PF  ;
  302.  
  303. *
  304. **** Les flux aux interfaces sont dans de
  305. *    CHPOINT FACE
  306. *
  307.  
  308. FLUX1  = 'EXCO' TFLUX 'RN';
  309. FLUX2X = 'EXCO' TFLUX 'RUX';
  310. FLUX2Y = 'EXCO' TFLUX 'RUY';
  311. FLUX3  = 'EXCO' TFLUX 'RETN';
  312.  
  313. FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE)) ;
  314. FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE)) ;
  315.  
  316. *
  317. **** Test 'KONV'
  318. *
  319.  
  320. f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
  321. ERRO = 'ABS' (1D-12 '*' f1gd);
  322. LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO;
  323.  
  324. f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
  325. ERRO = 'ABS' (1D-12 '*' f2gd);
  326. LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO;
  327. LOGI1 = LOGI1 'ET' LOGI2;
  328.  
  329. f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
  330. ERRO = 'ABS' (1D-11 '*' f3gd) ;
  331. LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO;
  332. LOGI1 = LOGI1 'ET' LOGI2;
  333.  
  334. f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
  335. ERRO = 'ABS' (1D-12 '*' f4gd) ;
  336. LOGI2 = ('ABS' ((f4 '*' ORIENT) '-' f4gd)) < ERRO;
  337. LOGI1 = LOGI1 'ET' LOGI2;
  338.  
  339. 'SI' ('NON' LOGI1);
  340.     'MESSAGE' ;
  341.     'MESSAGE' 'OPERATEUR KONV';
  342.     'MESSAGE' ('CHAINE' METO);
  343.     'MESSAGE' ;
  344.     'MESSAGE' ('CHAINE' 'df1 ' ('ABS' (f1gd '-'  (ORIENT '*'f1)))
  345.                         ' f1 ' f1);                       
  346.     'MESSAGE' ('CHAINE' 'df2 ' ('ABS' (f2gd '-'  (ORIENT '*'f2)))                       
  347.                         ' f2 ' f2);                       
  348.     'MESSAGE' ('CHAINE' 'df3 ' ('ABS' (f3gd '-'  (ORIENT '*'f3)))                       
  349.                         ' f3 ' f3);                       
  350.     'MESSAGE' ('CHAINE' 'df4 ' ('ABS' (f4gd '-'  (ORIENT '*'f4)))                           
  351.                         ' f4 ' f4);
  352.     'ERREUR' 5;
  353. 'FINSI' ;
  354.  
  355. ************************************************
  356. ***** Splitting des scalaires passifs **********
  357. ************************************************
  358.  
  359. TFLUX2 DT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' 'VLH'
  360.   $DOMTOT PGAZ2 ('MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'TUTU' 'TATA')
  361.   ROF VITF PF SCAF ;
  362.  
  363. ERRO = ( 'MAXIMUM' (TFLUX '-'
  364.          ('EXCO'  ('MOTS' 'RN' 'RUX' 'RUY' 'RETN') TFLUX2)) 'ABS') ;
  365.  
  366. 'SI' (ERRO > 1.0D-15);
  367.     'MESSAGE' 'Probleme dans les transport de scalaires passifs' ;
  368.     'ERREUR' 5;
  369. 'FINSI' ;
  370.  
  371. *
  372. **** Test 'KONV' 
  373. *
  374.  
  375. fceltutu = 'EXTRAIRE' TFLUX2  'TUTU' P1;
  376. LOGI1 = ('ABS' ((fceltutu '*' ORIENT) '-' ftutugd)) <
  377.         ('ABS' (1.0D-12 * fceltutu)) ;
  378.  
  379. fceltata = 'EXTRAIRE' TFLUX2 'TATA' P1;
  380. LOGI2 = ('ABS' ((fceltata '*' ORIENT) '-' ftatagd)) <
  381.         ('ABS' (1.0D-12 * fceltata)) ;
  382. LOGI1 = LOGI1 'ET' LOGI2;
  383.  
  384. 'SI' ('NON' LOGI1);
  385.     'MESSAGE' ;
  386.     'MESSAGE' 'OPERATEUR KONV';
  387.     'MESSAGE' 'Transport de scalaires passifs' ;
  388.     'ERREUR' 5;
  389. 'FINSI' ;
  390.  
  391. **************************
  392. **** 'KONV', TEST SS  ****
  393. **************************
  394.  
  395. TFLUX DT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' 'SS'
  396.   $DOMTOT PGAZ ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')  ROF VITF PF  ;
  397.  
  398. *
  399. **** Les flux aux interfaces sont dans de
  400. *    CHPOINT FACE
  401. *
  402.  
  403. FLUX1  = 'EXCO' TFLUX 'RN';
  404. FLUX2X = 'EXCO' TFLUX 'RUX';
  405. FLUX2Y = 'EXCO' TFLUX 'RUY';
  406. FLUX3  = 'EXCO' TFLUX 'RETN';
  407.  
  408. FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE)) ;
  409. FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE)) ;
  410.  
  411. *
  412. **** Test 'KONV'
  413. *
  414.  
  415. f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
  416. ERRO = 'ABS' (1D-12 '*' f1gd);
  417. LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO;
  418.  
  419. f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
  420. ERRO = 'ABS' (1D-12 '*' f2gd);
  421. LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO;
  422. LOGI1 = LOGI1 'ET' LOGI2;
  423.  
  424. f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
  425. ERRO = 'ABS' (1D-11 '*' f3gd) ;
  426. LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO;
  427. LOGI1 = LOGI1 'ET' LOGI2;
  428.  
  429. f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
  430. ERRO = 'ABS' (1D-12 '*' f4gd) ;
  431. LOGI2 = ('ABS' ((f4 '*' ORIENT) '-' f4gd)) < ERRO;
  432. LOGI1 = LOGI1 'ET' LOGI2;
  433.  
  434. 'SI' ('NON' LOGI1);
  435.     'MESSAGE' ;
  436.     'MESSAGE' 'OPERATEUR KONV';
  437.     'MESSAGE' ('CHAINE' METO);
  438.     'MESSAGE' ;
  439.     'MESSAGE' ('CHAINE' 'df1 ' ('ABS' (f1gd '-'  (ORIENT '*'f1)))
  440.                         ' f1 ' f1);                       
  441.     'MESSAGE' ('CHAINE' 'df2 ' ('ABS' (f2gd '-'  (ORIENT '*'f2)))                       
  442.                         ' f2 ' f2);                       
  443.     'MESSAGE' ('CHAINE' 'df3 ' ('ABS' (f3gd '-'  (ORIENT '*'f3)))                       
  444.                         ' f3 ' f3);                       
  445.     'MESSAGE' ('CHAINE' 'df4 ' ('ABS' (f4gd '-'  (ORIENT '*'f4)))                           
  446.                         ' f4 ' f4);
  447.     'ERREUR' 5;
  448. 'FINSI' ;
  449.  
  450. ************************************************
  451. ***** Splitting des scalaires passifs **********
  452. ************************************************
  453.  
  454. TFLUX2 DT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' 'SS'
  455.   $DOMTOT PGAZ2 ('MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'TUTU' 'TATA')
  456.   ROF VITF PF SCAF ;
  457.  
  458. ERRO = ( 'MAXIMUM' (TFLUX '-'
  459.          ('EXCO'  ('MOTS' 'RN' 'RUX' 'RUY' 'RETN') TFLUX2)) 'ABS') ;
  460.  
  461. 'SI' (ERRO > 1.0D-15);
  462.     'MESSAGE' 'Probleme dans les transport de scalaires passifs' ;
  463.     'ERREUR' 5;
  464. 'FINSI' ;
  465.  
  466. *
  467. **** Test 'KONV' 
  468. *
  469.  
  470. fceltutu = 'EXTRAIRE' TFLUX2  'TUTU' P1;
  471. LOGI1 = ('ABS' ((fceltutu '*' ORIENT) '-' ftutugd)) <
  472.         ('ABS' (1.0D-12 * fceltutu)) ;
  473.  
  474. fceltata = 'EXTRAIRE' TFLUX2 'TATA' P1;
  475. LOGI2 = ('ABS' ((fceltata '*' ORIENT) '-' ftatagd)) <
  476.         ('ABS' (1.0D-12 * fceltata)) ;
  477. LOGI1 = LOGI1 'ET' LOGI2;
  478.  
  479. 'SI' ('NON' LOGI1);
  480.     'MESSAGE' ;
  481.     'MESSAGE' 'OPERATEUR KONV';
  482.     'MESSAGE' 'Transport de scalaires passifs' ;
  483.     'ERREUR' 5;
  484. 'FINSI' ;
  485.  
  486. ****************************************************
  487. ****************************************************
  488. ******** Fin boucle sur les angles         *********
  489. ****************************************************
  490. ****************************************************
  491.  
  492. 'FIN' BLOC;
  493.  
  494. 'FIN' ;
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  
  507.  

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