Accueil | Cast3M

* fichier :  konv_ther_sup.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 en cas de supersonicite               ****
****                                                   ****
**** Methodes: VLH                                     ****
****                                                   ****
**** A. BECCANTINI DRN/DMT/SEMT/TTMF    DECEMBRE 1998  ****
***********************************************************
 
'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' ;
 
*
**** Ordre des polynoms cv_i
*
 
PGAZ . 'NORD' = 4 ;
 
*
**** Especes qui sont dans les equations d'Euler (rien)
*
 
*
**** Espece qui n'y est pas
*
 
 
PGAZ . 'ESPNEULE' = 'H2  ';
 
*
 
PGAZ .  'H2  ' = 'TABLE'  ;
 
*
**** R (J/Kg/K)
*
 
PGAZ .  'H2  ' . 'R' = 4130.0 ;
 
*
**** Regressions polynomials
*
 
PGAZ .  '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 ;
 
 
*************************** 
***** 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 
*
 
rog = 1.412 ;
pg = 1.17D5 ;
ung  = 1180.34 ;
utg  = 18.2 ;
*  cg = 350
 
rod = 1.0 ;
pd = 4.1117D5 ;
und  = 1188.34 ;
utd  = 1.2 ;
 
*
**** tg
*
 
tg = pg '/' (rog * (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 = tg ;
T2 = T * T ;
T3 = T2 * T ;
T4 = T3 * T ;
T5 = T4 * T ;
 
 
 ETG = (A0H2 * T) '+' (A1H2 * T2 '/' 2) '+' (A2H2 * T3 '/' 3) '+'
      (A3H2 * T4 '/' 4) '+' (A4H2 * T5 '/' 5) ;
 
 ecing = 0.5D0 '*'  ((ung '*' ung) '+' (utg '*' utg));
 
 retg = rog '*' (ETG '+' ecing) ;
 
*
***  flux en (n,t)
*
 
f1gd = ung  '*'  rog;
f2gd = (rog '*' (ung '*' ung)) '+' pg ;
f3gd = rog '*' ung '*' utg ;
f4gd = ung '*' (retg '+' pg);
 
 
****************************************************
****************************************************
******** Boucle sur les angles             *********
****************************************************
****************************************************
 
DANGLE = 360 '/' 7.15;
ANGLE = 11.3 ;
 
'REPETER' BLMETO 2 ;
 
'SI' ('EGA' &BLMETO 1) ;
  METO = 'VLH'  ;
'FINSI' ;
 
'SI' ('EGA' &BLMETO 2) ;
 METO =  'SS'  ;
'FINSI' ;
 
'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;
 
    'ELIMINATION' (DOM1 ET DOM2) 1D-6;
    DOMTOT = DOM1 ET DOM2;
 
 $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  ****
***********************
 
uxg = (ung '*' ('COS' ANGLE)) '-' (utg '*' ('SIN' ANGLE));
uyg = (ung '*' ('SIN' ANGLE)) '+' (utg '*' ('COS' ANGLE));
uxd = (und '*' ('COS' ANGLE)) '-' (utd '*' ('SIN' ANGLE));
uyd = (und '*' ('SIN' ANGLE)) '+' (utd '*' ('COS' ANGLE));
 
 
RN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' rog;
RN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' rod;
RN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
 
VN1   = 'KCHT' $DOM1   'VECT' 'CENTRE' (uxg uyg);
VN2   = 'KCHT' $DOM2   'VECT' 'CENTRE' (uxd uyd);
VN    = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VN1 'ET' VN2);
 
PN1   = 'KCHT' $DOM1   'SCAL' 'CENTRE' pg;
PN2   = 'KCHT' $DOM2   'SCAL' 'CENTRE' pd;
PN    = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
 
 
***************************
****  L'operateur  PRET****
***************************
 
ORDESP = 1;
ORDTEM = 1;
 
 
ROF VITF PF  =  'PRET' 'PERFTEMP'  ORDESP ORDTEM
    $DOMTOT PGAZ RN VN PN ;
 
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 = 'SIGN' ORIENT;   
 
* 'KONV'
 
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' ;
 
TFLUX DELTAT =  'KONV' 'VF' 'PERFTEMP' 'FLUX' METO
         $DOMTOT PGAZ LISTINCO  ROF VITF PF ;
 
*
**** Les flux aux interfaces sont dans de
*    CHPOINT FACE
*
 
FLUX1 = 'EXCO' 'RN' TFLUX ;
FLUX2X = 'EXCO' 'RUX' TFLUX ;
FLUX2Y = 'EXCO' 'RUY' TFLUX ;
FLUX3 = 'EXCO' 'RETN' TFLUX ;
 
FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));
 
 
*
**** Test 'KONV'
*
 
f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
ERRO = 'ABS' (1D-12 '*' f1gd);
LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO;
 
f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
ERRO = 'ABS' (1D-12 '*' f2gd);
LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
 
f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
ERRO = 'ABS' (1D-11 '*' f3gd) ;
LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
 
f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
ERRO = 'ABS' (1D-12 '*' 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' ;
 
 
****************************************************
****************************************************
******** Fin boucle sur les angles         *********
****************************************************
****************************************************
 
'FIN' BLOC;
 
'FIN' BLMETO ;
 
'FIN' ;