* fichier : tube2D.dgibi
******************************************************************
* CALCUL DU TUBE A CHOC 2D *
* FORMULATION VF COMPRESSIBLE EXPLICITE *
* DIFFERENTS SOLVEURS *
* *
* A. BECCANTINI TTMF MARS 1998 *
******************************************************************
'OPTION' 'DIME' 2 ;
'OPTION' 'ELEM' 'QUA4' ;
'OPTION' 'ISOV' 'SULI' ;
'OPTION' 'ECHO' 1 ;
'OPTION' 'TRAC' 'X' ;
GRAPH = VRAI ;
GRAPH = FAUX ;
************
* MAILLAGE *
************
RAF = 4 ;
NX=25 '*' RAF ;
NY= RAF ;
L1 = 1. ;
H1 = (L1 '/' NX) * NY ;
A3 = 0.0 0.0 ;
A4 = 0.0 H1 ;
A3A4 = A3 'DROIT' NY A4 ;
A1A2 = A3A4 'PLUS' ((-1 * L1) 0.0) ;
A5A6 = A3A4 'PLUS' (L1 0.0) ;
DOM1 = A1A2 'REGLER' NX A3A4 ;
DOM2 = A3A4 'REGLER' NX A5A6 ;
DOMTOT = DOM1 'ET' DOM2 ;
$DOMTOT = 'MODE' DOMTOT 'EULER' ;
$DOM1 = 'MODE' DOM1 'EULER' ;
$DOM2 = 'MODE' DOM2 'EULER' ;
TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
TDOM1 = 'DOMA' $DOM1 'VF' ;
TDOM2 = 'DOMA' $DOM2 'VF' ;
MDOMTOT = TDOMTOT . 'QUAF' ;
MDOM1 = TDOM1 . 'QUAF' ;
MDOM2 = TDOM2 . 'QUAF' ;
'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) ((H1 '/' NY) '/' 100.) ;
*
******* Creation de la ligne Utilisée pour le Post-Traitement
* reliant les points centres
*
PINI = (TDOMTOT . 'CENTRE') 'POIN' 'PROC' ((-1 * L1) 0.0) ;
PFIN = (TDOMTOT . 'CENTRE') 'POIN' 'PROC' (L1 0.0) ;
IAUX = (2 * NX) - 1 ;
COURB = PINI 'DROIT' IAUX PFIN;
'ELIMINATION' Courb ('DOMA' $DOMTOT 'CENTRE')
((H1 '/' NY) '/' 100.) ;
************************************************************
* CONDITIONS INITIALES ET LIMITES. *
************************************************************
gamgd = 1.4D0;
*
*** Etat gauche
*
rog = 1.0 ;
pg = 1.0 ;
reg = pg '/' (gamgd '-' 1.0) ;
*
*** Etat droite
*
rod = 1.0D-1 ;
pd = 1.0D-1;
red = (pd '/' (gamgd '-' 1.0)) ;
*
*
*** gamma
*
GAMN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gamgd;
*
*** ro
*
RO_1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog;
RO_2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RO_1 'ET' RO_2) ;
*
*** ro u, ro v
*
GN_1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (0.0 0.0);
GN_2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (0.0 0.0);
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN_1 'ET' GN_2);
*
*** ro e
*
RE_1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' reg ;
RE_2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' red ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RE_1 ET RE_2) ;
********************************************************
*** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM ***
********************************************************
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE';
*
**** Les debits dans le repaire tube
*
GNX = 'EXCO' 'UX' GN;
GNY = 'EXCO' 'UY' GN;
*
*** GRAPHIQUE DES C.I.
*
'SI' GRAPH ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ;
CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' 0.0);
TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' 0.0);
TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'GN at t=' 0.0);
'FINSI' ;
*
**** Les variables primitives
*
VN PN = 'PRIM' 'PERFMONO'
RN GN RETN GAMN ;
*
**** Les gradients
*
TITI CHLIM COEFG = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN ;
*
*************************************************************************
*************************************************************************
*********************** LE CALCUL ***************************************
*************************************************************************
*************************************************************************
*
*** Methodes possibles :
*
* 'VANLEER'
* 'VLH '
* 'HUSVLH '
* 'GODUNOV'
* 'AUSMPLUS'
* 'HLLC '
* ...
METO = 'HLLC' ;
LOGSO = VRAI ;
LOGRK2 = VRAI ;
SAFFAC = 0.7 ;
NITER = 1000 ;
TFINAL = 0.2 ;
*
TPS = 0.0 ;
ZERO = 1.0D-12 ;
'TEMPS' 'ZERO' ;
'REPETER' BL1 NITER ;
*
**** Primitive variables
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
'SI' LOGSO ;
GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN 'GRADGEO' COEFG ;
GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN 'GRADGEO' COEFG ;
GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN 'GRADGEO' COEFG ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
$DOMTOT
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
GAMN ;
'SINON' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
$DOMTOT
RN
VN
PN
GAMN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')
ROF VITF PF GAMF ;
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
*
**** Total time step
*
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
*
**** Increment of the variables (convection)
*
SI LOGRK2 ;
RESIDU = 0.5 * DTMIN '*' RESIDU ;
DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU
('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
*
RNP = RN '+' DRN ;
GNP = GN '+' DGN ;
RETNP = RETN '+' DRETN ;
VNP PNP = 'PRIM' 'PERFMONO' RNP GNP RETNP GAMN ;
'SI' LOGSO ;
GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RNP 'GRADGEO' COEFG ;
GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PNP 'GRADGEO' COEFG ;
GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VNP 'GRADGEO' COEFG ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
$DOMTOT
RNP GRADR ALR
VNP GRADV ALV
PNP GRADP ALP
GAMN ;
'SINON' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
$DOMTOT
RNP
VNP
PNP
GAMN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')
ROF VITF PF GAMF ;
'FINSI' ;
RESIDU = DTMIN '*' RESIDU ;
DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU
('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
TPS = TPS '+' DTMIN ;
RN = RN '+' DRN ;
GN = GN '+' DGN ;
RETN = RETN '+' DRETN ;
'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ;
'FINSI' ;
'SI' (TPS > TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
'TEMPS' ;
*
*************************************************************************
*************************************************************************
*********************** FIN CALCUL *************************************
*************************************************************************
*************************************************************************
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
*
**** La vitesse dans le repaire tube
*
VNN = 'EXCO' 'UX' VN;
VNT = 'EXCO' 'UY' VN;
*
*** GRAPHIQUE DES SOLUTIONS
*
'SI' (GRAPH);
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ;
CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' TPS);
TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' TPS);
TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'RUN at t=' TPS);
'FINSI' ;
*
*** Objects evolutions
*
'SI' LOGSO ;
IE = 2 ;
'SINON' ;
IE = 1 ;
'FINSI' ;
SI LOGRK2 ;
IT = 2 ;
'SINON' ;
IT = 1 ;
'FINSI' ;
xx yy = 'COORDONNEE' Courb;
ss = xx ;
evxx = 'EVOL' 'CHPO' ss Courb ;
lxx = 'EXTRAIRE' evxx 'ORDO' ;
x0 = 'MINIMUM' lxx;
x1 = 'MAXIMUM' lxx;
evro = 'EVOL' 'CHPO' RN Courb ;
lro = 'EXTRAIRE' evro 'ORDO' ;
evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro;
tit = ('CHAINE' '1D ' METO ' IT ' IT
' IE ' IE ' TPS ' TPS) ;
evu = 'EVOL' 'CHPO' VNN Courb ;
lu = 'EXTRAIRE' evu 'ORDO' ;
evu = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
evv = 'EVOL' 'CHPO' VNT Courb ;
lv = 'EXTRAIRE' evv 'ORDO' ;
evv = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
SI GRAPH;
'DESSIN' evv 'TITRE' tit 'XBOR' x0 x1;
FINSI;
evp = 'EVOL' 'CHPO' PN Courb ;
lp = 'EXTRAIRE' evp 'ORDO' ;
evp = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
ls = lro '**' gamgd;
ls = lp '/' ls;
evs = 'EVOL' 'MANU' 'x' lxx 's' ls ;
ts = CHAINE '1D ' METO ' : s IT ' IT ' IE ' IE
' tmps ' TPS ' elem ' 'CUB8' ;
*
*** Solution analytique
*
lxxa = 'PROG' -0.49000 -0.47000 -0.45000 -0.43000 -0.41000 -0.390
-0.37000 -0.35000 -0.33000 -0.31000 -0.29000 -0.27000 -0.25000
-0.23000 -0.21000 -0.19000 -0.17000 -0.15000 -0.13000 -0.11000
-9.00000E-2 -7.00000E-2 -5.00000E-2 -3.00000E-2 -1.00000E-2;
lxxa = lxxa 'ET' ('PROG'
1.00000E-02 3.00000E-02 5.00000E-02 7.00000E-02 9.00000E-02
0.11000 0.13000 0.15000 0.17000 0.19000 0.21000 0.23000 0.25000
0.27000 0.29000 0.31000 0.33000 0.35000 0.37000 0.39000 0.41000
0.43000 0.45000 0.47000) ;
lpa = 'PROG' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99193
.85405 .79151 .69598 .64297 .57841 .51997 .45356 .41754
.37307 .33114 .29475 .28482 .28482 .28482 .28482 .28482
.28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482
.28482 .28482 .28482 .28482 .28482 .28482 .10000 .10000
.10000 .10000 .10000;
lroa = 'PROG' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99423
.89343 .84619 .77192 .72945 .67635 .62680 .56851 .53588
.49447 .45410 .41786 .40776 .40776 .40776 .40776 .40776
.40776 .40776 .40776 .40776 .40776 .20444 .20444 .20444
.20444 .20444 .20444 .20444 .20444 .20444 .10000 .10000
.10000 .10000 .10000;
lua = 'PROG' 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 6.84663E-03 .13185 .19435 .29851 .36173 .44507 .52768
.63185 .69395 .77728 .86406 .94740 .97167 .97167 .97167
.97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167
.97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167
0.0 0.0 0.0 0.0 0.0;
EPSI=0.26;
lsa = lroa '**' gamgd;
lsa = lpa '/' lsa;
evroa = 'EVOL' 'MANU' 'x' lxxa 'ro' lroa;
evpa = 'EVOL' 'MANU' 'x' lxxa 'p' lpa;
evua = 'EVOL' 'MANU' 'x' lxxa 'ua' lua;
evsa = 'EVOL' 'MANU' 'x' lxxa 's' lsa;
*PM du fait de la précision du SIN et du COS, (de l'ordre de 10-10)
* l'intervalle lxx est plus petit que lxxa.
* Comme IPOL fait une erreur dans ce cas, on restreint lxxa
'REMPLACER' lxxa 1 ('EXTR' lxx 1) ;
'REMPLACER' lxxa ('DIME' lxxa) ('EXTR' lxx ('DIME' lxx)) ;
lper = 'IPOL' lxxa lxx lp ;
luer = 'IPOL' lxxa lxx lu ;
lroer = 'IPOL' lxxa lxx lro ;
dlp = 'MAXIMUM' (ABS ( lper '-' lpa));
dlu = 'MAXIMUM' (ABS ( luer '-' lua));
dlro = 'MAXIMUM' (ABS ( lroer '-' lroa));
dl = (dlp '+' dlu '+' dlro) '/' 3.0D0;
*
*** Quelque DESSIN
*
TAB1=TABLE;
TAB1.'TITRE'= TABLE ;
TAB1.1='MARQ TRIB REGU';
TAB1.'TITRE' . 1 = MOT 'Numerique';
TAB1.2='MARQ CROI REGU';
TAB1.'TITRE' . 2 = MOT 'Exacte';
SI GRAPH;
'DESSIN' (evro 'ET' evroa) 'XBOR' x0 x1
'LEGE' TAB1 'TITRE' tit ;
'DESSIN' (evp 'ET' evpa) 'XBOR' x0 x1
'LEGE' TAB1 'TITRE' tit ;
'DESSIN' (evu 'ET' evua) 'XBOR' x0 x1
'LEGE' TAB1 'TITRE' tit ;
'DESSIN' (evs 'ET' evsa) 'XBOR' x0 x1
'LEGE' TAB1 'TITRE' tit ;
FINSI;
SI (DL > EPSI);
MESSAGE 'Erreur calcul du tube a choc';
MESSAGE dl;
ERREUR 5;
FINSI;
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer