* fichier : tube3Daxi.dgibi
******************************************************************
* CALCUL DU TUBE A CHOC *
* FORMULATION VF COMPRESSIBLE EXPLICITE *
* DIFFERENTS SOLVEURS *
* 3DAXIS *
* *
* A. BECCANTINI LTMF FEVRIER 2004 *
******************************************************************
'OPTION' 'DIME' 3 'ELEM' 'CUB8' 'ISOV' 'SULI'
'ECHO' 1 'TRAC' 'X' ;
GRAPH = VRAI ;
GRAPH = FAUX ;
*
*** Methodes possibles :
*
* 'VANLEER'
* 'VLH '
* 'HUSVLH '
* 'GODUNOV'
* 'AUSMPLUS'
METO = 'VLH' ;
LOGSO = VRAI ;
SAFFAC = 0.5 ;
NITER = 100000 ;
TFINAL = 0.2 ;
************
* MAILLAGE *
************
RAF = 1 ;
NX= 50 '*' RAF ;
L = 1.0D0 ;
DX = 0.5 '*' L '/' NX ;
AM1 = 0.0 (-0.5 '*' L) 0.0 ;
A0 = 0.0 0.0 0.0 ;
A1 = 0.0 (0.5 '*' L) 0.0 ;
BM1 = DX (-0.5 '*' L) 0.0 ;
B0 = DX 0.0 0.0 ;
B1 = DX (0.5 '*' L) 0.0 ;
*
A0B0 = A0 'DROIT' 1 B0 ;
A1B1 = A1 'DROIT' 1 B1 ;
AM1BM1 = AM1 'DROIT' 1 BM1 ;
*
BAS1 = A0B0 'REGLER' NX AM1BM1 ;
BAS2 = A0B0 'REGLER' NX A1B1 ;
*
BAS1M = BAS1 'TOURNER' -5. (0.0 0.0 0.0) (0.0 1.0 0.0) ;
BAS1P = BAS1 'TOURNER' 5. (0.0 0.0 0.0) (0.0 1.0 0.0) ;
BAS2M = BAS2 'TOURNER' -5. (0.0 0.0 0.0) (0.0 1.0 0.0) ;
BAS2P = BAS2 'TOURNER' 5. (0.0 0.0 0.0) (0.0 1.0 0.0) ;
'ELIMINATION' (BAS1M 'ET' BAS2M) (DX '/' 1000.) ;
'ELIMINATION' (BAS1P 'ET' BAS2P) (DX '/' 1000.) ;
DOM1 = 'VOLUME' BAS1M 1 BAS1P ;
DOM2 = 'VOLUME' BAS2M 1 BAS2P ;
'ELIMINATION' DOM1 (DX '/' 1000.) ;
DOM1 = 'REGENERER' DOM1 ;
'ELIMINATION' DOM2 (DX '/' 1000.) ;
DOM2 = 'REGENERER' DOM2 ;
DOMTOT = DOM1 'ET' DOM2 ;
'ELIMINATION' DOMTOT (DX '/' 1000.) ;
$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) (DX '/' 100.) ;
*
******* Creation de la ligne Utilisée pour le Post-Traitement
* reliant les points centres
*
CENCOR = 'COORDONNEE' 2 (TDOMTOT . 'CENTRE') ;
PINI = 'POIN' (CENCOR 'POIN' 'MINI') 1 ;
PFIN = 'POIN' (CENCOR 'POIN' 'MAXI') 1 ;
IAUX = (2 * NX) - 1 ;
COURB = PINI 'DROIT' IAUX PFIN;
'ELIMINATION' (1.0D-6) Courb ('DOMA' $DOMTOT 'CENTRE');
'SI' GRAPH ;
'TRACER' (DOMTOT 'ET' (COURB 'COULEUR' ROUG))
'TITR' '3D domain' ;
'FINSI' ;
************************************************************
* CONDITIONS INITIALES ET LIMITES. *
************************************************************
gamgd = 1.4D0;
*
*** Etat gauche
*
rog = 1.0 ;
pg = 1.0 ;
retg = (pg '/' (gamgd '-' 1.0)) ;
*
*** Etat droite
*
rod = 1.0D-1 ;
pd = 1.0D-1;
retd = (pd '/' (gamgd '-' 1.0)) ;
*
*** gamma
*
GAMN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gamgd;
*
*** ro
*
RN = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' rog
'NATU' 'DISCRET') 'ET' ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1
'SCAL' rod 'NATU' 'DISCRET') ;
*
*** ro u, ro v
*
GN = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 3 'UX' 0.0 'UY' 0.0
'UZ' 0.0 'NATU' 'DISCRET') ;
*
*** ro e
*
RETN = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' retg
'NATU' 'DISCRET') 'ET' ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1
'SCAL' retd 'NATU' 'DISCRET') ;
********************************************************
*** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM ***
********************************************************
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE';
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
*
'SI' GRAPH ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ;
'TRACER' CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' 0.0);
'TRACER' CHM_PN MOD1 'TITR' ('CHAINE' 'P at t=' 0.0);
'TRACER' CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' 0.0);
'FINSI' ;
*
**** Les gradients
*
CACCA CHLIM COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN ;
CACCA CHLIM COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
('MOTS' 'UX' 'UY' 'UZ') VN ;
*
**** Le calcul
*
TPS = 0.0 ;
ZERO = 1.0D-12 ;
DTMIN = 0.0 ;
'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' COEFSCAL ;
GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN 'GRADGEO' COEFSCAL ;
GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY' 'UZ') VN 'GRADGEO' COEFVECT ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
$DOMTOT
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
GAMN ;
* ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 2
* $DOMTOT
* RN GRADR ALR
* VN GRADV ALV
* PN GRADP ALP
* GAMN (DTMIN '/' 2) ;
'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' 'RUZ' '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)
*
RESIDU = DTMIN '*' RESIDU ;
DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY' 'RUZ') RESIDU
('MOTS' 'UX' 'UY' 'UZ') ;
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 ' DTPS =' DTMIN
' TPS =' TPS) ;
'FINSI' ;
'SI' (TPS '>EG' TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
'TEMPS' ;
*
***** On calcule les variables primitive
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
*
**** La vitesse dans le repaire tube
*
VNX = 'EXCO' 'UX' VN;
VNY = 'EXCO' 'UY' VN;
VNZ = 'EXCO' 'UZ' VN;
VNN = VNY ;
VNT = VNX ;
VNV = VNZ ;
GNN = VNN * RN;
GNT = VNT * RN;
GNV = VNV * RN;
*
*** GRAPHIQUE DES SOLUTIONS
*
'SI' (GRAPH);
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN;
CHM_GNN = 'KCHA' $DOMTOT 'CHAM' GNN ;
CHM_GNT = 'KCHA' $DOMTOT 'CHAM' GNT ;
CHM_GNV = 'KCHA' $DOMTOT 'CHAM' GNV ;
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_GNN MOD1 'TITR' ('CHAINE' 'ROUN at t=' TPS);
TRAC CHM_GNT MOD1 'TITR' ('CHAINE' 'ROUT at t=' TPS);
'FINSI' ;
*
*** Objects evolutions
*
'SI' LOGSO ;
IE = 2 ;
'SINON' ;
IE = 1 ;
'FINSI' ;
IT = 1 ;
SS = 'COORDONNEE' 2 Courb;
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;
tro = CHAINE '1D ' METO ' : RO IT ' IT ' IE ' IE
' tmps ' TPS ' elem ' 'CUB8' ;
evu = 'EVOL' 'CHPO' VNN Courb ;
lu = 'EXTRAIRE' evu 'ORDO' ;
evu = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
tu = CHAINE '1D ' METO ' : U IT ' IT ' IE ' IE
' tmps ' TPS ' elem ' 'CUB8' ;
evv = 'EVOL' 'CHPO' VNT Courb ;
lv = 'EXTRAIRE' evv 'ORDO' ;
evv = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
tv = CHAINE '1D ' METO ' : V IT ' IT ' IE ' IE
' tmps ' TPS ' elem ' 'CUB8' ;
SI GRAPH;
'DESSIN' evv 'TITRE' tv 'XBOR' x0 x1;
FINSI;
evw = 'EVOL' 'CHPO' VNV Courb ;
lw = 'EXTRAIRE' evw 'ORDO' ;
evw = 'EVOL' 'MANU' 'x' lxx 'w' lw ;
tw = CHAINE '1D ' METO ' : W IT ' IT ' IE ' IE
' tmps ' TPS ' elem ' 'CUB8' ;
SI GRAPH;
'DESSIN' evw 'TITRE' tw 'XBOR' x0 x1;
FINSI;
evp = 'EVOL' 'CHPO' PN Courb ;
lp = 'EXTRAIRE' evp 'ORDO' ;
evp = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
tp = CHAINE '1D ' METO ' : P IT ' IT ' IE ' IE
' tmps ' TPS ' elem ' 'CUB8' ;
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;
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) 'TITRE' tro 'XBOR' x0 x1
'LEGE' TAB1;
'DESSIN' (evp 'ET' evpa) 'TITRE' tp 'XBOR' x0 x1
'LEGE' TAB1;
'DESSIN' (evu 'ET' evua) 'TITRE' tu 'XBOR' x0 x1
'LEGE' TAB1;
'DESSIN' (evs 'ET' evsa) 'TITRE' ts 'XBOR' x0 x1
'LEGE' TAB1;
FINSI;
SI (DLP > 1.0D-1);
MESSAGE 'Erreur calcul du tube a choc';
MESSAGE dl;
ERREUR 5;
FINSI;
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer