* fichier : konv_impl_murs.dgibi
***********************************************************
***********************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la ****
**** solution des ****
**** Equations d'Euler pour un gaz parfait ****
**** OPERATEURS PRIM, PRET, KONV ****
**** Implicit: calcul du jacobien du residu ****
**** ****
**** Cas gaz monoespece, "calorically perfect" ****
**** ****
**** Methodes: VLH ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/LTMF AOUT 2000 ****
***********************************************************
***********************************************************
'OPTION' 'DIME' 2 ;
'OPTION' 'ELEM' QUA4 ;
'OPTION' 'ECHO' 0 ;
'OPTION' 'TRAC' 'X' ;
*
*** GRAPH
*
GRAPH = FAUX ;
* GRAPH = VRAI ;
ERRTOL = 1.0D-3 ;
***************************
***** DOMAINE SPATIAL ****
***************************
A1 = 0.0D0 0.0D0;
A2 = 1.0D0 0.0D0;
A3 = 1.0D0 1.0D0;
A4 = 0.0D0 1.0D0;
L12 = A1 'DROIT' 1 A2;
L23 = A2 'DROIT' 1 A3;
L34 = A3 'DROIT' 1 A4;
L41 = A4 'DROIT' 1 A1;
DOM10 = 'DALL' L12 L23 L34 L41
'PLANE';
*
*** Etat
*
ro = 1.11 ;
p = 1234.7;
gam = 1.4D0;
cson = (gam '*' p '/' ro) '**' 0.5 ;
un = 0.2 * cson ;
ut = 0.3 * cson ;
*
*** retgd
*
ecin = 0.5D0 '*' ro '*' ((un '*' un ) '+' (ut '*' ut ));
ret = (p '/' (gam '-' 1.0)) '+' ecin ;
*
*** flux en (n,t)
*
****************************************************
****************************************************
******** Boucle sur les angles *********
****************************************************
****************************************************
DANGLE = (360.0 '/' 7.15) ;
ANGLE = 0.0 ;
'REPETER' BLOC 8 ;
*
*** Rotation
*
ANGLE = ANGLE '+' DANGLE ;
ORIG = 0.0D0 0.0D0;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Angle de rotation= ' ANGLE);
'MESSAGE' ;
DOMTOT = DOM10 'TOURNER' ANGLE ORIG;
$DOMTOT = 'MODELISER' DOMTOT 'EULER';
TDOMTOT = 'DOMA' $DOMTOT 'VF';
MDOMTOT = TDOMTOT . 'QUAF' ;
'ELIM' MDOMTOT 1.E-6 ;
PCON = ('DOMA' $DOMTOT 'CENTRE') 'POIN' 1 ;
'SI' GRAPH;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE' ) 'ET'
('DOMA' $DOMTOT 'CENTRE')) 'TITRE' 'Domaine et centre' ;
'FINSI' ;
***********************
**** Les CHPOINTs ****
***********************
ux = (un '*' ('COS' ANGLE)) '-' (ut '*' ('SIN' ANGLE));
uy = (un '*' ('SIN' ANGLE)) '+' (ut '*' ('COS' ANGLE));
ro0 = ro ;
gnx0 = ro '*' ux ;
gny0 = ro '*' uy ;
ret0 = ret ;
GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro0 ;
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx0 gny0) ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret0 ;
***************************
**** L'operateur PRIM****
***************************
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ;
***************************
**** L'operateur PRET****
***************************
ORDESP = 1;
ORDTEM = 1;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMAN ;
***************************
**** L'operateur KONV****
***************************
METO = 'VLH' ;
LISTINCO = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ;
CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
IJACO = 'KONV' 'VF' 'PERFMONO' 'JACOCONS'
$DOMTOT LISTINCO
METO RN VITESSE PRES GAMMAN ;
***********************
***** Le residu *******
***********************
DEBRN0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ;
DEBGNX0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ;
DEBGNY0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ;
DEBRETN0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ;
*
***** Le jacobien du residu et le residu
*
UNRO = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 1.0
'RUXN' 0.0 'RUYN' 0.0 'RETN' 0.0
'NATURE' 'DISCRET' ;
UNGX = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 0.0
'RUXN' 1.0 'RUYN' 0.0 'RETN' 0.0
'NATURE' 'DISCRET' ;
UNGY = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 0.0
'RUXN' 0.0 'RUYN' 1.0 'RETN' 0.0
'NATURE' 'DISCRET' ;
UNRET = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 0.0
'RUXN' 0.0 'RUYN' 0.0 'RETN' 1.0
'NATURE' 'DISCRET' ;
DDRHO = 'KOPS' IJACO 'MULT' UNRO ;
*
*
***** DDRHO contient: dRES_RN '/' dRN ('RN') ;
* dRES_GXN '/' dRN ('RUXN') ;
* dRES_GYN '/' dRN ('RUYN') ;
* dRES_RETN '/' dRN ('RETN') ;
* DDGX contient ...
*
*
DDGX = 'KOPS' IJACO 'MULT' UNGX ;
DDGY = 'KOPS' IJACO 'MULT' UNGY ;
DDRET = 'KOPS' IJACO 'MULT' UNRET ;
DRR = 'EXTRAIRE' DDRHO PCON 'RN' ;
DGXR = 'EXTRAIRE' DDRHO PCON 'RUXN' ;
DGYR = 'EXTRAIRE' DDRHO PCON 'RUYN' ;
DRETR = 'EXTRAIRE' DDRHO PCON 'RETN' ;
DRGX = 'EXTRAIRE' DDGX PCON 'RN' ;
DGXGX = 'EXTRAIRE' DDGX PCON 'RUXN' ;
DGYGX = 'EXTRAIRE' DDGX PCON 'RUYN' ;
DRETGX = 'EXTRAIRE' DDGX PCON 'RETN' ;
DRGY = 'EXTRAIRE' DDGY PCON 'RN' ;
DGXGY = 'EXTRAIRE' DDGY PCON 'RUXN' ;
DGYGY = 'EXTRAIRE' DDGY PCON 'RUYN' ;
DRETGY = 'EXTRAIRE' DDGY PCON 'RETN' ;
DRRET = 'EXTRAIRE' DDRET PCON 'RN' ;
DGXRET = 'EXTRAIRE' DDRET PCON 'RUXN' ;
DGYRET = 'EXTRAIRE' DDRET PCON 'RUYN' ;
DRETRET = 'EXTRAIRE' DDRET PCON 'RETN' ;
*********************************************************************
*********************************************************************
***** On calcule les residues pour ro1 = ro0 * (1.'+' DELTA) *******
*********************************************************************
*********************************************************************
DELTA = 1.0D-4 ;
ro1 = ro0 '*' (1.0 '+' DELTA) ;
gnx1 = gnx0 ;
gny1 = gny0 ;
ret1 = ret0 ;
GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ;
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ;
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMAN ;
CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
***********************
***** Le residu *******
***********************
DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ;
DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ;
DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ;
DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ;
*
**** On calcule le jacobien numeriquement
*
DRRN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ro0) ;
DGXRN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ro0) ;
DGYRN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ro0) ;
DRETRN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ro0) ;
ERR1 = ('ABS' (DRRN '-' DRR)) '*' (ro0 '/' (ro0 '*' cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 1';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGXRN '-' DGXR)) '*'
(ro0 '/' (ro0 '*' cson '*' cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 2';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGYRN '-' DGYR)) '*'
(ro0 '/' (ro0 '*' cson '*' cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 3';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DRETRN '-' DRETR)) '*'
(ro0 '/' (ret0 '*' cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 4';
'ERREUR' 5 ;
'FINSI' ;
***************************************************************************
***************************************************************************
***** On calcule les residues pour gnx1 = gnx0 '+' (DELTA ro0 cson) ******
***************************************************************************
***************************************************************************
DELTA = 1.0D-4 ;
ro1 = ro0 ;
gnx1 = gnx0 '+' (DELTA '*' ro0 '*' cson) ;
gny1 = gny0 ;
ret1 = ret0 ;
GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ;
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ;
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMAN ;
CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
***********************
***** Le residu *******
***********************
DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ;
DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ;
DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ;
DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ;
*
**** On calcule le jacobien numeriquement
*
DRGXN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ro0 '*' cson) ;
DGXGXN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ro0 '*' cson) ;
DGYGXN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ro0 '*' cson) ;
DRETGXN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ro0 '*' cson) ;
ERR1 = ('ABS' (DRGXN '-' DRGX)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 5';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGXGXN '-' DGXGX)) '/' cson ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 6';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGYGXN '-' DGYGX)) '/' cson ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 7';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DRETGXN '-' DRETGX)) '*'
(ro0 / ret0) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 8';
'ERREUR' 5 ;
'FINSI' ;
***************************************************************************
***************************************************************************
***** On calcule les residues pour gny1 = gny0 '+' (DELTA ro0 cson) ******
***************************************************************************
***************************************************************************
DELTA = 1.0D-4 ;
ro1 = ro0 ;
gnx1 = gnx0 ;
gny1 = gny0 '+' (DELTA '*' ro0 '*' cson) ;
ret1 = ret0 ;
GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ;
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ;
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMAN ;
CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
***********************
***** Le residu *******
***********************
DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ;
DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ;
DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ;
DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ;
*
**** On calcule le jacobien numeriquement
*
DRGYN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ro0 '*' cson) ;
DGXGYN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ro0 '*' cson) ;
DGYGYN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ro0 '*' cson) ;
DRETGYN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ro0 '*' cson) ;
ERR1 = ('ABS' (DRGYN '-' DRGY)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 9';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGXGYN '-' DGXGY)) '/' cson ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 10';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGYGYN '-' DGYGY)) '/' cson ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 11';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DRETGYN '-' DRETGY)) '*'
(ro0 / ret0) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 12';
'ERREUR' 5 ;
'FINSI' ;
***************************************************************************
***************************************************************************
***** On calcule les residues pour ret1 = ret0 '*' (1. '+' DELTA) ******
***************************************************************************
***************************************************************************
DELTA = 1.0D-4 ;
ro1 = ro0 ;
gnx1 = gnx0 ;
gny1 = gny0 ;
ret1 = ret0 * (1. '+' DELTA) ;
GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ;
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ;
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMAN ;
CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
***********************
***** Le residu *******
***********************
DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ;
DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ;
DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ;
DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ;
*
**** On calcule le jacobien numeriquement
*
DRRETN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ret0) ;
DGXRETN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ret0) ;
DGYRETN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ret0) ;
DRETRETN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ret0) ;
ERR1 = ('ABS' (DRRETN '-' DRRET)) '*' (ret0 '/' (ro0 '*' cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 13';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGXRETN '-' DGXRET)) '*'
(ret0 '/' (ro0 '*' cson * cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 14';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DGYRETN '-' DGYRET)) '*'
(ret0 '/' (ro0 '*' cson * cson)) ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 15';
'ERREUR' 5 ;
'FINSI' ;
ERR1 = ('ABS' (DRETRETN '-' DRETRET)) '/' cson ;
'SI' (ERR1 > ERRTOL) ;
'MESSAGE' 'Problem 16';
'ERREUR' 5 ;
'FINSI' ;
'FIN' BLOC ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer