Download konv_resi_ther_cons.dgibi

Back to the list

* fichier :  konv_resi_ther_cons.dgibi
***********************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la    ****
**** solution des                                      **** 
**** Equations d'Euler pour un gaz parfait             ****
**** OPERATEURS PRET, KONV                             ****
****                                                   ****
**** Cas gaz monoespece  "thermally perfect"           ****
**** Consistence                                       ****
****                                                   ****
**** Methodes: VLH/SS                                  ****
****                                                   ****
**** A. BECCANTINI DM2S/SEMT/LTMF    AVRIL 2001        ****
***********************************************************


 'OPTION'  'DIME' 2 ;
 'OPTION'  'ELEM' QUA4 ;
 'OPTION'  'ECHO'  0 ;
 'OPTION'  'TRAC' 'X' ;

*
*** GRAPH
*

 GRAPH = FAUX ;
* GRAPH = VRAI ;

**** Cas monoespece : la table proprieté de gaz

 PGAZ = 'TABLE' ;
 PGAZ2 = 'TABLE' ;

*
**** Ordre des polynoms cv_i
*

 PGAZ . 'NORD' = 4 ;
 PGAZ2 . 'NORD' = 4 ;

*
**** Especes qui sont dans les equations d'Euler (rien)
*

*
**** Espece qui n'y est pas
*


 PGAZ . 'ESPNEULE' = 'H2  ';
 PGAZ2 . 'ESPNEULE' = 'H2  ';

*

 PGAZ .  'H2  ' = 'TABLE'  ;
 PGAZ2 .  'H2  ' = 'TABLE'  ;

*
**** R (J/Kg/K)
*

 PGAZ .  'H2  ' . 'R' = 4130.0 ;
 PGAZ2 .  'H2  ' . 'R' = 4130.0 ;

*
**** Regressions polynomials
* 

 PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
                               -2.37281455E-07 1.84701105E-11 ;
 PGAZ2 .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
                               -2.37281455E-07 1.84701105E-11 ;
                              
*
**** "Enthalpies" (ou energies) de formations a OK (J/Kg)
*     Note: ce sont des entites numeriques
*

 PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
 PGAZ2 .  'H2  ' . 'H0K' = -4.195D6 ;


* PGAZ2: transport de scalaires passifs

 PGAZ2 . 'SCALPASS' = 'MOTS' 'TUTU' 'TATA' ;

*************************** 
***** DOMAINE SPATIAL  ****
***************************


 A1 = 0.0D0 0.0D0;
 A2 = 1.0D0 0.0D0;
 A3 = 2.0D0 0.0D0;
 A4 = 2.0D0 1.0D0;
 A5 = 1.0D0 1.0D0;
 A6 = 0.0D0 1.0D0;

 L12 = A1 'DROIT' 1 A2;
 L23 = A2 'DROIT' 1 A3;
 L34 = A3 'DROIT' 1 A4;
 L45 = A4 'DROIT' 1 A5;
 L56 = A5 'DROIT' 1 A6;
 L61 = A6 'DROIT' 1 A1;
 L25 = A2 'DROIT' 1 A5;


 DOM10  = 'DALL' L12 L25 L56 L61
        'PLANE';
 DOM20  = 'DALL' L23 L34 L45 ('INVERSE' L25)
        'PLANE';
*
*** Point ou on controlle la consistence
*

 P10 = 1.0 0.5;

*
*** Etats gauche et droite 
*

 rogd = 1.412 ;
 pgd = 1.17D6 ;
 ungd  = 8.34 ;
 utgd  = 18.2 ;
 tutugd = 14.9 ;
 tatagd = 23.8 ;

*
**** tg
*

 tgd = pgd '/' (rogd * (PGAZ .  'H2  ' . 'R')) ;


*
**** gamg, gamd , htg , htd
*
*

 A0H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 1 ;
 A1H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 2 ;
 A2H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 3 ;
 A3H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 4 ;
 A4H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 5 ;

 T = tgd ;
 T2 = T * T ;
 T3 = T2 * T ;
 T4 = T3 * T ;
 T5 = T4 * T ;

 cvgd =  A0H2 '+' (A1H2 * T) '+' (A2H2 * T2) '+' (A3H2 * T3) '+'
        (A4H2 * T4) ;
        
 ETGd = (A0H2 * T) '+' (A1H2 * T2 '/' 2) '+' (A2H2 * T3 '/' 3) '+'
      (A3H2 * T4 '/' 4) '+' (A4H2 * T5 '/' 5) ;

 ecingd = 0.5D0 '*'  ((ungd '*' ungd) '+' (utgd '*' utgd));
           
 retgd = rogd '*' (ETGd '+' ecingd) ;

 gamgd = (cvgd '+' (PGAZ .  'H2  ' . 'R')) '/' cvgd ;
  
*
***  flux en (n,t)
*

 f1gd = ungd  '*'  rogd ;
 f2gd = (rogd '*' (ungd '*' ungd)) '+' pgd ;
 f3gd = rogd '*' ungd '*' utgd ;
 f4gd = ungd '*' (retgd '+' pgd);
 ftutugd =  ungd  '*'  rogd '*' tutugd ;
 ftatagd =  ungd  '*'  rogd '*' tatagd ;



****************************************************
****************************************************
******** Boucle sur les angles             *********
****************************************************
****************************************************

 DANGLE = 360 '/' 7.15;
 ANGLE = 11.3 ;

'REPETER' BLOC 8;

*
*** Rotation
*


    ANGLE = ANGLE '+' DANGLE;
    ORIG = 0.0D0 0.0D0;

    'MESSAGE' ;
    'MESSAGE'  (CHAIN 'Angle de rotation= ' ANGLE);
    'MESSAGE' ;

    DOM1   = DOM10   'TOURNER' ANGLE ORIG;
    DOM2   = DOM20   'TOURNER' ANGLE ORIG;
    P1     = P10     'TOURNER' ANGLE ORIG;

    DOMTOT = DOM1 ET DOM2;
    'ELIMINATION' DOMTOT 1D-6;

    $DOMTOT = 'MODELISER' DOMTOT 'EULER';
    $DOM1 = 'MODELISER'  DOM1 'EULER';
    $DOM2 = 'MODELISER'  DOM2 'EULER';
    TDOMTOT = 'DOMA'  $DOMTOT 'VF';
    TDOM1 = 'DOMA'  $DOM1 'VF';
    TDOM2 = 'DOMA'  $DOM2 'VF';
    MDOM1 = TDOM1 . 'QUAF' ;
    MDOM2 = TDOM2 . 'QUAF' ;
    MDOMTOT = TDOMTOT . 'QUAF' ;

    'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) 1.E-6 ;

    'SI' GRAPH;
       'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET'  
       ('DOMA' $DOMTOT 'FACEL') 'ET' P1) 'TITRE' 'Domaine et FACEL';
    'FINSI' ;

*
**** Redefinition de P1 dans $DOMTOT 'FACE'
*

    P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;

***********************
**** Les CHPOINTs  ****
***********************

    uxgd = (ungd '*' ('COS' ANGLE)) '-' (utgd '*' ('SIN' ANGLE));
    uygd = (ungd '*' ('SIN' ANGLE)) '+' (utgd '*' ('COS' ANGLE));

    RN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' rogd;
    RN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' rogd;;
    RN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);

    VN1   = 'KCHT' $DOM1   'VECT' 'CENTRE' (uxgd uygd);
    VN2   = 'KCHT' $DOM2   'VECT' 'CENTRE' (uxgd uygd);
    VN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VN1 'ET' VN2);

    PN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' pgd;
    PN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' pgd;
    PN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);

    SCAN1   = 'KCHT' $DOM1   'VECT' 'CENTRE'
              'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ;
    SCAN2   = 'KCHT' $DOM2   'VECT' 'CENTRE'
              'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ;
    SCAN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
              'COMP' (PGAZ2 . 'SCALPASS') (SCAN1 'ET' SCAN2);

***************************
****  L'operateur  PRET****
***************************

    ORDESP = 1;
    ORDTEM = 1;


    ROF VITF PF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
       $DOMTOT PGAZ RN VN PN ;

    ROF2 VITF2 PF2 SCAF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
       $DOMTOT PGAZ2 RN VN PN SCAN ;

    GEOP1    = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
    GEOP2    = ('DOMA' $DOMTOT 'FACE')  'ELEM' 'APPUYE' 'LARGEMENT' P1;


    REFGEOP1 = 'REDU' VITF GEOP2;

*
**** Orientation de la normal n de castem par raport a la
*    notre; t est par consequence
*

    NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1;
    NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1;

    'SI' (('ABS' NCOS) > ('ABS' NSIN));
       ORIENT = ('COS' ANGLE) '/' NCOS;
    'SINON'; 
       ORIENT = ('SIN' ANGLE) '/' NSIN;
    'FINSI' ;

    ORIENT = 'SIGNE' 'FLOTTANT' ORIENT;   



********************************
********************************
**** 'KONV'                 ****
********************************
********************************


    'REPETER' BLMETO 2 ;

       'SI' ('EGA' &BLMETO 1) ;
           METO = 'VLH'  ;
       'FINSI' ;

       'SI' ('EGA' &BLMETO 2) ;
           METO = 'SS'  ;
       'FINSI' ;

       'MESSAGE' ('CHAINE' 'METO = ' METO) ;
       LISTINC1 = ('MOTS' 'RN' 'RUNX' 'RUNY' 'RETN') ;
       CHPFLU DT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' METO
          $DOMTOT PGAZ LISTINC1 ROF VITF PF  ;

       FLUX1 = 'EXCO' ('EXTRAIRE' LISTINC1 1) CHPFLU ;
       FLUX2X = 'EXCO' ('EXTRAIRE' LISTINC1 2) CHPFLU ;
       FLUX2Y = 'EXCO' ('EXTRAIRE' LISTINC1 3) CHPFLU ;
       FLUX3 = 'EXCO' ('EXTRAIRE' LISTINC1 4) CHPFLU ;

       FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
       FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));


       f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
       ERRO = 1D-8 '*' f1gd;
       LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO ; 

       f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
       ERRO = 1D-8 '*' f2gd;
       LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO;
       LOGI1 = LOGI1 'ET' LOGI2;

       f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
       ERRO = 1D-8 '*' f3gd;
       LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO;
       LOGI1 = LOGI1 'ET' LOGI2;

       f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
       ERRO = 1D-8 '*' f4gd;
       LOGI2 = ('ABS' ((f4 '*' ORIENT) '-' f4gd)) < ERRO;
       LOGI1 = LOGI1 'ET' LOGI2;

       'SI' ('NON' LOGI1);
          'MESSAGE' ;
          'MESSAGE' 'OPERATEUR KONV';
          'MESSAGE' ('CHAINE' METO);
          'MESSAGE' ;
          'MESSAGE' ('CHAINE' 'df1 ' ('ABS' (f1gd '-'  (ORIENT '*'f1)))
                     ' f1 ' f1);                       
          'MESSAGE' ('CHAINE' 'df2 ' ('ABS' (f2gd '-'  (ORIENT '*'f2)))                       
                        ' f2 ' f2);                       
          'MESSAGE' ('CHAINE' 'df3 ' ('ABS' (f3gd '-'  (ORIENT '*'f3)))                       
                        ' f3 ' f3);                       
          'MESSAGE' ('CHAINE' 'df4 ' ('ABS' (f4gd '-'  (ORIENT '*'f4)))                           
                        ' f4 ' f4);
          'ERREUR' 5;
       'FINSI' ;

*
**** Scalaires passifs
*                 

       'MESSAGE' ('CHAINE' 'METO = ' METO ' (scalaires passifs)') ;
       LISTINC2 = ('MOTS' 'RN' 'RUNX' 'RUNY' 'RETN' 'TUTU' 'TATA') ;
       CHPFLU2 DT2 = 'KONV' 'VF' 'PERFTEMP' 'FLUX' METO
          $DOMTOT PGAZ2 LISTINC2 ROF VITF PF SCAF ;

       CHPFLU3 = 'EXCO' LISTINC1 CHPFLU2 LISTINC1 ;

       'SI' (('MAXIMUM' (CHPFLU3 '-' CHPFLU) 'ABS') > 1.0D-15) ;
          'MESSAGE' 'Probleme dans les transport de scalaires passifs' ;
          'ERREUR' 5;
       'FINSI' ;

       fceltutu = 'EXTRAIRE' CHPFLU2  'TUTU' P1;
       LOGI1 = ('ABS' ((fceltutu '*' ORIENT) '-' ftutugd)) <
               ('ABS' (1.0D-8 * fceltutu)) ;
        
       fceltata = 'EXTRAIRE' CHPFLU2 'TATA' P1;
       LOGI2 = ('ABS' ((fceltata '*' ORIENT) '-' ftatagd)) <
                 ('ABS' (1.0D-8 * fceltata)) ;
       LOGI1 = LOGI1 'ET' LOGI2;

       'SI' ('NON' LOGI1);
          'MESSAGE' ;
          'MESSAGE' 'OPERATEUR KONV';
          'MESSAGE' 'Transport de scalaires passifs' ;
          'ERREUR' 5;
       'FINSI' ;


*
****** Le residu
*
*
****** CHPDIR = chpoint sur les face de chaque elts
*
*
*      NB Dans ce cas particulier, DOMTOT n'a que de QUA4,
*         donc  ($DOMTOT . 'XXNORMAE') contient une seule
*         zona elementaire!

      
       RESI MATCEL = 'KOPS' 'MATRIK' ;
       MAICHA = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAE') 'MAILLAGE' ;
       NBELEM = 'NBEL' MAICHA ;
       
       'REPETER' BLELE NBELEM ;
          CHPDIR MATCEL = 'KOPS' 'MATRIK' ;
          ELEM1 = &BLELE 'ELEM' MAICHA ;
          ELEMC1 = 'CHANGER' ELEM1 'POI1' ;
          NBPOIN = 'NBNO' ELEMC1 ;
          PCEN = ('DOMA' $DOMTOT 'CENTRE') 'POIN' &BLELE ;
          XVOL = ('DOMA' $DOMTOT 'XXVOLUM') 'EXTRAIRE' PCEN 'SCAL' ;
          MAICEN = 'MANU' 'POI1'PCEN ;
          'REPETER' BLEPO NBPOIN ;
             P1 = &BLEPO 'POIN' ELEMC1 ;
             ELEP1 = 'MANUEL' 'POI1' P1 ;
             COEFNO = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAE')
                'SCAL' 1 &BLELE &BLEPO ;
             CHPDIR = CHPDIR 'ET' ('MANUEL' 'CHPO' ELEP1 1 'SCAL'
                COEFNO 'NATU' 'DISCRET') ;
          'FIN' BLEPO ;
          MAIDIR = 'EXTRAIRE' CHPDIR 'MAILLAGE' ;
*
********* Le flux qui sort
*          
          FLUSOR = CHPDIR '*' ('REDU' CHPFLU2 MAIDIR)
             ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL' 'SCAL' 'SCAL')
              LISTINC2 LISTINC2 ;
          CHPUN = 'MANUEL' 'CHPO' MAIDIR 1 'SCAL' (-1./XVOL) ;
          NBCOMP = 'DIME' LISTINC2 ;
          'REPETER' BLCOM NBCOMP ;
              NOMCOM = 'EXTRAIRE' LISTINC2 &BLCOM ;
              CELSCA = 'XTY' FLUSOR CHPUN ('MOTS' NOMCOM)
                 ('MOTS' 'SCAL') ;
              RESI = RESI 'ET' ('MANUEL' 'CHPO' MAICEN 1
              NOMCOM CELSCA 'NATU' 'DISCRET') ;
          'FIN' BLCOM ;
       'FIN' BLELE ;   
                      
       CHPRES DT1 = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
          $DOMTOT PGAZ2 LISTINC2 ROF VITF PF SCAF ;
                
       ERRO = 'MAXIMUM' (RESI '-' CHPRES) 'ABS' ;
       'SI' (ERRO > 1.0D-8) ;
          'MESSAGE' ;
          'MESSAGE' 'Probleme dans KONV' ;
          'ERREUR' 5 ;
       'FINSI' ;

       ERRO = 'ABS' (DT '-' DT1) ;
       'SI' (ERRO > 1.0D-15) ;
          'MESSAGE' ;
          'MESSAGE' 'Probleme dans KONV' ;
          'ERREUR' 5 ;
       'FINSI' ;

       

   'FIN' BLMETO ;       

****************************************************
****************************************************
******** Fin boucle sur les angles         *********
****************************************************
****************************************************

'FIN' BLOC;

'FIN' ;

 

 

 

© Cast3M 2003 - All rights reserved.
Disclaimer