Télécharger la plaquette Cast3M


fr en

Télécharger konv_resi_ther_cons.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  konv_resi_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. **** Consistence                                       ****
  12. ****                                                   ****
  13. **** Methodes: VLH/SS                                  ****
  14. ****                                                   ****
  15. **** A. BECCANTINI DM2S/SEMT/LTMF    AVRIL 2001        ****
  16. ***********************************************************
  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.     DOMTOT = DOM1 ET DOM2;
  208.     'ELIMINATION' DOMTOT 1D-6;
  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.     RN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' rogd;
  241.     RN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' rogd;;
  242.     RN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
  243.  
  244.     VN1   = 'KCHT' $DOM1   'VECT' 'CENTRE' (uxgd uygd);
  245.     VN2   = 'KCHT' $DOM2   'VECT' 'CENTRE' (uxgd uygd);
  246.     VN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VN1 'ET' VN2);
  247.  
  248.     PN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' pgd;
  249.     PN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' pgd;
  250.     PN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
  251.  
  252.     SCAN1   = 'KCHT' $DOM1   'VECT' 'CENTRE'
  253.               'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ;
  254.     SCAN2   = 'KCHT' $DOM2   'VECT' 'CENTRE'
  255.               'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ;
  256.     SCAN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
  257.               'COMP' (PGAZ2 . 'SCALPASS') (SCAN1 'ET' SCAN2);
  258.  
  259. ***************************
  260. ****  L'operateur  PRET****
  261. ***************************
  262.  
  263.     ORDESP = 1;
  264.     ORDTEM = 1;
  265.  
  266.  
  267.     ROF VITF PF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
  268.        $DOMTOT PGAZ RN VN PN ;
  269.  
  270.     ROF2 VITF2 PF2 SCAF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
  271.        $DOMTOT PGAZ2 RN VN PN SCAN ;
  272.  
  273.     GEOP1    = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
  274.     GEOP2    = ('DOMA' $DOMTOT 'FACE')  'ELEM' 'APPUYE' 'LARGEMENT' P1;
  275.  
  276.  
  277.     REFGEOP1 = 'REDU' VITF GEOP2;
  278.  
  279. *
  280. **** Orientation de la normal n de castem par raport a la
  281. *    notre; t est par consequence
  282. *
  283.  
  284.     NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1;
  285.     NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1;
  286.  
  287.     'SI' (('ABS' NCOS) > ('ABS' NSIN));
  288.        ORIENT = ('COS' ANGLE) '/' NCOS;
  289.     'SINON'; 
  290.        ORIENT = ('SIN' ANGLE) '/' NSIN;
  291.     'FINSI' ;
  292.  
  293.     ORIENT = 'SIGN' ORIENT;   
  294.  
  295.  
  296.  
  297. ********************************
  298. ********************************
  299. **** 'KONV'                 ****
  300. ********************************
  301. ********************************
  302.  
  303.  
  304.     'REPETER' BLMETO 2 ;
  305.  
  306.        'SI' ('EGA' &BLMETO 1) ;
  307.            METO = 'VLH'  ;
  308.        'FINSI' ;
  309.  
  310.        'SI' ('EGA' &BLMETO 2) ;
  311.            METO = 'SS'  ;
  312.        'FINSI' ;
  313.  
  314.        'MESSAGE' ('CHAINE' 'METO = ' METO) ;
  315.        LISTINC1 = ('MOTS' 'RN' 'RUNX' 'RUNY' 'RETN') ;
  316.        CHPFLU DT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' METO
  317.           $DOMTOT PGAZ LISTINC1 ROF VITF PF  ;
  318.  
  319.        FLUX1 = 'EXCO' ('EXTRAIRE' LISTINC1 1) CHPFLU ;
  320.        FLUX2X = 'EXCO' ('EXTRAIRE' LISTINC1 2) CHPFLU ;
  321.        FLUX2Y = 'EXCO' ('EXTRAIRE' LISTINC1 3) CHPFLU ;
  322.        FLUX3 = 'EXCO' ('EXTRAIRE' LISTINC1 4) CHPFLU ;
  323.  
  324.        FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
  325.        FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));
  326.  
  327.  
  328.        f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
  329.        ERRO = 1D-8 '*' f1gd;
  330.        LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO ; 
  331.  
  332.        f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
  333.        ERRO = 1D-8 '*' f2gd;
  334.        LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO;
  335.        LOGI1 = LOGI1 'ET' LOGI2;
  336.  
  337.        f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
  338.        ERRO = 1D-8 '*' f3gd;
  339.        LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO;
  340.        LOGI1 = LOGI1 'ET' LOGI2;
  341.  
  342.        f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
  343.        ERRO = 1D-8 '*' f4gd;
  344.        LOGI2 = ('ABS' ((f4 '*' ORIENT) '-' f4gd)) < ERRO;
  345.        LOGI1 = LOGI1 'ET' LOGI2;
  346.  
  347.        'SI' ('NON' LOGI1);
  348.           'MESSAGE' ;
  349.           'MESSAGE' 'OPERATEUR KONV';
  350.           'MESSAGE' ('CHAINE' METO);
  351.           'MESSAGE' ;
  352.           'MESSAGE' ('CHAINE' 'df1 ' ('ABS' (f1gd '-'  (ORIENT '*'f1)))
  353.                      ' f1 ' f1);                       
  354.           'MESSAGE' ('CHAINE' 'df2 ' ('ABS' (f2gd '-'  (ORIENT '*'f2)))                       
  355.                         ' f2 ' f2);                       
  356.           'MESSAGE' ('CHAINE' 'df3 ' ('ABS' (f3gd '-'  (ORIENT '*'f3)))                       
  357.                         ' f3 ' f3);                       
  358.           'MESSAGE' ('CHAINE' 'df4 ' ('ABS' (f4gd '-'  (ORIENT '*'f4)))                           
  359.                         ' f4 ' f4);
  360.           'ERREUR' 5;
  361.        'FINSI' ;
  362.  
  363. *
  364. **** Scalaires passifs
  365. *                 
  366.  
  367.        'MESSAGE' ('CHAINE' 'METO = ' METO ' (scalaires passifs)') ;
  368.        LISTINC2 = ('MOTS' 'RN' 'RUNX' 'RUNY' 'RETN' 'TUTU' 'TATA') ;
  369.        CHPFLU2 DT2 = 'KONV' 'VF' 'PERFTEMP' 'FLUX' METO
  370.           $DOMTOT PGAZ2 LISTINC2 ROF VITF PF SCAF ;
  371.  
  372.        CHPFLU3 = 'EXCO' LISTINC1 CHPFLU2 LISTINC1 ;
  373.  
  374.        'SI' (('MAXIMUM' (CHPFLU3 '-' CHPFLU) 'ABS') > 1.0D-15) ;
  375.           'MESSAGE' 'Probleme dans les transport de scalaires passifs' ;
  376.           'ERREUR' 5;
  377.        'FINSI' ;
  378.  
  379.        fceltutu = 'EXTRAIRE' CHPFLU2  'TUTU' P1;
  380.        LOGI1 = ('ABS' ((fceltutu '*' ORIENT) '-' ftutugd)) <
  381.                ('ABS' (1.0D-8 * fceltutu)) ;
  382.  
  383.        fceltata = 'EXTRAIRE' CHPFLU2 'TATA' P1;
  384.        LOGI2 = ('ABS' ((fceltata '*' ORIENT) '-' ftatagd)) <
  385.                  ('ABS' (1.0D-8 * fceltata)) ;
  386.        LOGI1 = LOGI1 'ET' LOGI2;
  387.  
  388.        'SI' ('NON' LOGI1);
  389.           'MESSAGE' ;
  390.           'MESSAGE' 'OPERATEUR KONV';
  391.           'MESSAGE' 'Transport de scalaires passifs' ;
  392.           'ERREUR' 5;
  393.        'FINSI' ;
  394.  
  395.  
  396. *
  397. ****** Le residu
  398. *
  399. *
  400. ****** CHPDIR = chpoint sur les face de chaque elts
  401. *
  402. *
  403. *      NB Dans ce cas particulier, DOMTOT n'a que de QUA4,
  404. *         donc  ($DOMTOT . 'XXNORMAE') contient une seule
  405. *         zona elementaire!
  406.  
  407.  
  408.        RESI MATCEL = 'KOPS' 'MATRIK' ;
  409.        MAICHA = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAE') 'MAILLAGE' ;
  410.        NBELEM = 'NBEL' MAICHA ;
  411.  
  412.        'REPETER' BLELE NBELEM ;
  413.           CHPDIR MATCEL = 'KOPS' 'MATRIK' ;
  414.           ELEM1 = &BLELE 'ELEM' MAICHA ;
  415.           ELEMC1 = 'CHANGER' ELEM1 'POI1' ;
  416.           NBPOIN = 'NBNO' ELEMC1 ;
  417.           PCEN = ('DOMA' $DOMTOT 'CENTRE') 'POIN' &BLELE ;
  418.           XVOL = ('DOMA' $DOMTOT 'XXVOLUM') 'EXTRAIRE' PCEN 'SCAL' ;
  419.           MAICEN = 'MANU' 'POI1'PCEN ;
  420.           'REPETER' BLEPO NBPOIN ;
  421.              P1 = &BLEPO 'POIN' ELEMC1 ;
  422.              ELEP1 = 'MANUEL' 'POI1' P1 ;
  423.              COEFNO = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAE')
  424.                 'SCAL' 1 &BLELE &BLEPO ;
  425.              CHPDIR = CHPDIR 'ET' ('MANUEL' 'CHPO' ELEP1 1 'SCAL'
  426.                 COEFNO 'NATU' 'DISCRET') ;
  427.           'FIN' BLEPO ;
  428.           MAIDIR = 'EXTRAIRE' CHPDIR 'MAILLAGE' ;
  429. *
  430. ********* Le flux qui sort
  431. *          
  432.           FLUSOR = CHPDIR '*' ('REDU' CHPFLU2 MAIDIR)
  433.              ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL' 'SCAL' 'SCAL')
  434.               LISTINC2 LISTINC2 ;
  435.           CHPUN = 'MANUEL' 'CHPO' MAIDIR 1 'SCAL' (-1./XVOL) ;
  436.           NBCOMP = 'DIME' LISTINC2 ;
  437.           'REPETER' BLCOM NBCOMP ;
  438.               NOMCOM = 'EXTRAIRE' LISTINC2 &BLCOM ;
  439.               CELSCA = 'XTY' FLUSOR CHPUN ('MOTS' NOMCOM)
  440.                  ('MOTS' 'SCAL') ;
  441.               RESI = RESI 'ET' ('MANUEL' 'CHPO' MAICEN 1
  442.               NOMCOM CELSCA 'NATU' 'DISCRET') ;
  443.           'FIN' BLCOM ;
  444.        'FIN' BLELE ;   
  445.  
  446.        CHPRES DT1 = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
  447.           $DOMTOT PGAZ2 LISTINC2 ROF VITF PF SCAF ;
  448.  
  449.        ERRO = 'MAXIMUM' (RESI '-' CHPRES) 'ABS' ;
  450.        'SI' (ERRO > 1.0D-8) ;
  451.           'MESSAGE' ;
  452.           'MESSAGE' 'Probleme dans KONV' ;
  453.           'ERREUR' 5 ;
  454.        'FINSI' ;
  455.  
  456.        ERRO = 'ABS' (DT '-' DT1) ;
  457.        'SI' (ERRO > 1.0D-15) ;
  458.           'MESSAGE' ;
  459.           'MESSAGE' 'Probleme dans KONV' ;
  460.           'ERREUR' 5 ;
  461.        'FINSI' ;
  462.  
  463.  
  464.  
  465.    'FIN' BLMETO ;       
  466.  
  467. ****************************************************
  468. ****************************************************
  469. ******** Fin boucle sur les angles         *********
  470. ****************************************************
  471. ****************************************************
  472.  
  473. 'FIN' BLOC;
  474.  
  475. 'FIN' ;
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  

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