Télécharger la plaquette Cast3M


fr en

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 = 'SIGN' 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.  
  454.  
  455.  

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