* fichier : sinebum_fmm.dgibi
*********************************************************************
* *
* SINE-SHAPED BUMP *
* *
* CALCUL DE L ECOULEMENT SUBSONIQUE ISENTROPIQUE STATIONNAIRE DANS *
* UN CANAL *
* *
* Methode implicite sans matrice *
* Boundary conditions imposed via ghost cells *
* *
* BECCANTINI A., SFME/LTMF, MAI 2002 *
* *
*********************************************************************
TYEL = 'QUA4' ;
'OPTION' 'DIME' 2 'ELEM' TYEL 'ISOV' 'SULI'
'ECHO' 0 'TRAC' 'X' ;
GRAPH = VRAI ;
GRAPH = FAUX ;
******************
******************
******************
**** MAILLAGE ****
******************
******************
******************
******************
*
*
*
RAF = 2 ;
NY = 5 '*' RAF ;
NX1 = 4 '*' RAF ;
NX2 = 2 '*' NX1 ;
NX3 = NX1 ;
NX = (NX1 '+' NX2 '+' NX3) ;
DX = (4.0 '/' NX) ;
A0 = -2.0 0.0 ;
A1 = -1.0 0.0 ;
A2 = 1.0 0.0 ;
A3 = 2.0 0.0 ;
A4 = 2.0 1.0 ;
A5 = -2.0 1.0 ;
*
**** LIGB
*
LIGB1 = A0 'DROIT' NX1 A1 ;
* LIGB2 (On utilise un propriete de 'ET' ; si 'ET' change ?)
xcel = ('COORDONNEE' 1 A1) '+' DX ;
ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
ACEL = xcel ycel ;
LIGB2 = A1 'DROIT' 1 ACEL ;
'REPETER' BL1 (NX2 '-' 2) ;
ACEL0 = ACEL ;
xcel = xcel '+' DX ;
ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
ACEL = xcel ycel ;
LIGB2 = LIGB2 'ET' (ACEL0 'DROIT' 1 ACEL) ;
'FIN' BL1;
LIGB2 = LIGB2 'ET' (ACEL 'DROIT' 1 A2) ;
LIGB3 = A2 'DROIT' NX3 A3 ;
LIGB = LIGB1 'ET' LIGB2 'ET' LIGB3 ;
*
**** LIGH
*
LIGH = A4 'DROIT' NX A5 ;
*
**** DOMINT
*
DOMINT = LIGH 'REGLER' NY ('INVERSE' LIGB) ;
LIGCON = 'CONTOUR' DOMINT ;
*
*** LIGG
*
LIGG = LIGCON 'ELEM' 'COMP' A5 A0 ;
*
**** LIGD
*
LIGD = LIGCON 'ELEM' 'COMP' A3 A4 ;
FROD = LIGD 'TRANSLATION' 1 (DX 0.0) ;
FROG = LIGG 'TRANSLATION' 1 ((-1.0 '*' DX) 0.0) ;
* DOMTOT = total region
DOMTOT = DOMINT 'ET' FROD 'ET' FROG ;
'ELIMINATION' (1.0D-3 '/' RAF) DOMTOT ;
MDOMTOT = 'MODELISER' DOMTOT 'EULER' ;
MDOMINT = 'MODELISER' DOMINT 'EULER' ;
MFROD = 'MODELISER' FROD 'EULER' ;
MFROG = 'MODELISER' FROG 'EULER' ;
MLIGG = 'MODELISER' LIGG 'EULER' ;
MLIGD = 'MODELISER' LIGD 'EULER' ;
*
**** Creation of DOMAINE tables via the MODEL object
*
TDOMTOT = 'DOMA' MDOMTOT 'VF' ;
TDOMINT = 'DOMA' MDOMINT 'VF' ;
TFROD = 'DOMA' MFROD 'VF' ;
TFROG = 'DOMA' MFROG 'VF' ;
TLIGG = 'DOMA' MLIGG 'VF' ;
TLIGD = 'DOMA' MLIGD 'VF' ;
QDOMTOT = TDOMTOT . 'QUAF' ;
QDOMINT = TDOMINT . 'QUAF' ;
QFROD = TFROD . 'QUAF' 'COULEUR' 'ROUG' ;
QFROG = TFROG . 'QUAF' 'COULEUR' 'ROUG' ;
QLIGD = TLIGD . 'QUAF' ;
QLIGG = TLIGG . 'QUAF' ;
'ELIMINATION' QDOMTOT (1.0D-3 '/' RAF) QDOMINT ;
'ELIMINATION' QDOMTOT (1.0D-3 '/' RAF) QFROG ;
'ELIMINATION' QDOMTOT (1.0D-3 '/' RAF) QFROD ;
'ELIMINATION' QDOMTOT (1.0D-3 '/' RAF) QLIGG ;
'ELIMINATION' QDOMTOT (1.0D-3 '/' RAF) QLIGD ;
*
**** SEG2 mesh for BC
*
ELP1 = 'MANUEL' 'POI1' ('POIN' 1 ('DOMA' MLIGD 'CENTRE')) ;
ELTINT = ('DOMA' MDOMINT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT'
ELP1 ;
ELTFRO = 'ELEM' ('DOMA' MFROD 'FACEL') 'APPUYE' 'LARGEMENT'
ELP1 ;
P1 = 'POIN' 1 ELTINT ;
P2 = 'POIN' 1 ELTFRO ;
COND = 'MANUEL' 'SEG2' P1 P2 ;
NBL1 = ('NBNO' ('DOMA' MLIGD 'CENTRE')) '-' 1 ;
'REPETER' BL1 NBL1 ;
ELP1 = 'MANUEL' 'POI1' ('POIN' (&BL1 '+' 1)
('DOMA' MLIGD 'CENTRE')) ;
ELTINT = ('DOMA' MDOMINT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT'
ELP1 ;
ELTFRO = 'ELEM' ('DOMA' MFROD 'FACEL') 'APPUYE' 'LARGEMENT'
ELP1 ;
P1 = 'POIN' 1 ELTINT ;
P2 = 'POIN' 1 ELTFRO ;
COND = COND 'ET' ('MANUEL' 'SEG2' P1 P2) ;
'FIN' BL1 ;
COND = COND 'COULEUR' 'ROUG' ;
ELP1 = 'MANUEL' 'POI1' ('POIN' 1 ('DOMA' MLIGG 'CENTRE')) ;
ELTINT = ('DOMA' MDOMINT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT'
ELP1 ;
ELTFRO = 'ELEM' ('DOMA' MFROG 'FACEL') 'APPUYE' 'LARGEMENT'
ELP1 ;
P1 = 'POIN' 1 ELTINT ;
P2 = 'POIN' 1 ELTFRO ;
CONG = 'MANUEL' 'SEG2' P1 P2 ;
NBL1 = ('NBNO' ('DOMA' MLIGG 'CENTRE')) '-' 1 ;
'REPETER' BL1 NBL1 ;
ELP1 = 'MANUEL' 'POI1' ('POIN' (&BL1 '+' 1)
('DOMA' MLIGG 'CENTRE')) ;
ELTINT = ('DOMA' MDOMINT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT'
ELP1 ;
ELTFRO = 'ELEM' ('DOMA' MFROG 'FACEL') 'APPUYE' 'LARGEMENT'
ELP1 ;
P1 = 'POIN' 1 ELTINT ;
P2 = 'POIN' 1 ELTFRO ;
CONG = CONG 'ET' ('MANUEL' 'SEG2' P1 P2) ;
'FIN' BL1 ;
CONG = CONG 'COULEUR' 'ROUG' ;
'SI' FAUX ;
'TRACER' (DOMTOT 'ET' COND 'ET' CONG) 'TITRE'
'Domaine total' ;
'FINSI' ;
***************************************************************
***************************************************************
***************************************************************
************** Initial conditions *****************************
***************************************************************
***************************************************************
***************************************************************
GAMAIR = 1.4 ;
ro = 1.4 ;
p = 1.0 ;
u = 0.5;
***************************************************************
***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
***************************************************************
'DEBPROC' CONS ;
'ARGUMENT' RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
RVN = RN '*' VN ('MOTS' 'SCAL' 'SCAL') ('MOTS' 'UX' 'UY')
('MOTS' 'UX' 'UY') ;
CELL = 'PSCAL' RVN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
RECIN = 0.5 '*' CELL ;
REIN = PN '/' (GAMN '-' 1.0) ;
RETN = RECIN '+' REIN ;
DETR CELL ;
DETR RECIN ;
DETR REIN ;
'RESPRO' RVN RETN ;
'FINPROC' ;
***************************************************************
***************************************************************
***************************************************************
RN0 = 'MANU' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'SCAL' ro
'NATURE' 'DISCRET' ;
VN0 = ('MANU' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 2 'UX' u
'UY' 0.0 'NATURE' 'DISCRET') ;
PN0 = 'MANU' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'SCAL' p
'NATURE' 'DISCRET' ;
GAMN = 'MANU' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'SCAL' GAMAIR
'NATURE' 'DISCRET' ;
GN0 RETN0 = CONS RN0 VN0 PN0 GAMN ;
VN1 PN1 = 'PRIM' 'PERFMONO' RN0 GN0 RETN0 GAMN ;
ERRO = 'MAXIMUM' (PN1 '-' PN0) 'ABS' ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' 'Problem in the ic file!!!'
'ERREUR' 5 ;
'FINSI' ;
*
**** Plot of IC
*
MOD1 = 'MODELISER' ('DOMA' MDOMTOT 'MAILLAGE' ) 'THERMIQUE';
'SI' FAUX ;
CHM_RN = 'KCHA' MDOMTOT 'CHAM' RN0 ;
CHM_PN = 'KCHA' MDOMTOT 'CHAM' PN0 ;
CHM_VN = 'KCHA' MDOMTOT 'CHAM' VN0 ;
'TRACER' CHM_RN MOD1
'TITR' ('CHAINE' 'RN at t=' 0.0);
'TRACER' CHM_PN MOD1
'TITR' ('CHAINE' 'PN at t=' 0.0);
'TRACER' CHM_VN MOD1
'TITR' ('CHAINE' 'MACHN at t=' 0.0);
'FINSI' ;
MOTRN = 'RN' ;
MOTRNX = 'RUX' ;
MOTRNY = 'RUY' ;
MOTRETN = 'RETN' ;
LISTINCO = 'MOTS' MOTRN MOTRNX MOTRNY MOTRETN ;
**********************************************************
**********************************************************
**********************************************************
************** COMPUTATION OF THE SOLUTION ***************
**********************************************************
**********************************************************
**********************************************************
**********************************************************
ZERO = 1.0D-8 ;
*
**** Parameter for the computations
*
*
* Upwind scheme
* METO = 'VLH' ;
* METO = 'SS' ;
* METO = 'GODUNOV' ;
METO = 'AUSMPLUS' ;
* Reconstruction
* Iterations
* Final time
LOGREC = VRAI ;
NITER = 1000 ;
TFINAL = 10000.0D0 ;
* Overcoming of LCB (freezing of limiters)
NLCB = 1000 ;
TYPELIM = 'NOLIMITE' ;
* TYPELIM = 'LIMITEUR' ;
*
* Safety Factor fot the time step
* If SAFFAC < SAFFACM, then each FSAF-th iteration we compute
* SAFFAC = SAFFAC * 2 ;
*
SAFFAC = 10. ;
SAFFACM = 40. ;
FSAF = 400 ;
* Jacobi iterations
NJAC = 15 ;
EPSJAC = 1.0D-16 ;
RN = 'COPIER' RN0 ;
GN = 'COPIER' GN0 ;
RETN = 'COPIER' RETN0 ;
*
**** Coeff to compute gradients
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
GRADRN ALRN MCHRNC = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') RN ;
GRADPN ALPN MCHPNC = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') PN ;
GRADVN ALVN MCHVNC = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' TYPELIM
('MOTS' 'UX' 'UY') VN ;
* Limiters zero into ghost cells
LIMZER = ('MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' 1.0
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' MFROG 'CENTRE') 1 'SCAL' 0.0
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' MFROD 'CENTRE') 1 'SCAL' 0.0
'NATU' 'DISCRET') ;
* Names of the residuum components
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' ;
TPS = 0.0 ;
*
**** Temporal loop
*
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
LISTLINF = 'PROG' ;
LISTTPS = 'PROG' ;
LISTITER = 'PROG' ;
'TEMPS' ;
'REPETER' BL1 NITER ;
*
**** Primitive variables
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
'SI' LOGREC ;
GRADRN ALRN = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') RN 'GRADGEO' MCHRNC ;
GRADPN ALPN = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') PN 'GRADGEO' MCHPNC ;
GRADVN ALVN = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' TYPELIM
('MOTS' 'UX' 'UY') VN 'GRADGEO' MCHVNC ;
'SI' (&BL1 < NLCB) ;
* Limiters zero on ghost cells
ALRN0 = ALRN '*' LIMZER ;
ALPN0 = ALPN '*' LIMZER ;
ALVN0 = ALVN '*' LIMZER ;
'SINON' ;
'SI' (&BL1 'EGA' NLCB) ;
'MESSAGE' ;
'MESSAGE' 'On gele les limiteurs!!!' ;
'MESSAGE' ;
'FINSI' ;
'FINSI' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
MDOMTOT
RN GRADRN ALRN0
VN GRADVN ALVN0
PN GRADPN ALPN0
GAMN ;
'SINON' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
MDOMTOT
RN VN PN GAMN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
MDOMTOT LISTINCO ROF VITF PF GAMF ;
*
*** Time step
*
DTPS = DELTAT '*' SAFFAC ;
DTPS = 'MINIMUM' ('PROG' DTPS ((TFINAL '-' TPS) '*' 1.0001)) ;
TPS = TPS '+' DTPS ;
*
*** JACOBI
*
RNL0 = 'COPIER' RN ;
GNL0 = 'COPIER' GN ;
RETNL0 = 'COPIER' RETN ;
ANI = 'KONV' 'VF' 'PMONOFMM' 'AN' LISTINCO
MDOMTOT RN GN RETN GAMN SAFFAC ;
BNI = 'KONV' 'VF' 'PMONOFMM' 'BN' LISTINCO
MDOMTOT RN GN RETN GAMN ;
UNSANI = ('INVERSE' ANI) ;
DUN = UNSANI '*' RESIDU ;
'SI' VRAI ;
* On the ghost cells we impose DUN = 0
DUN = DUN '*' LIMZER ;
'FINSI' ;
DRN = 'EXCO' 'RN' DUN 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') DUN ('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' DUN 'SCAL' ;
RN = RNL0 '+' DRN ;
GN = GNL0 '+' DGN ;
RETN = RETNL0 '+' DRETN ;
'REPETER' BLJAC NJAC ;
DUN0 = DUN ;
BNP1 = 'KONV' 'VF' 'PMONOFMM' 'BN' LISTINCO
MDOMTOT RN GN RETN GAMN ;
RESP1 = RESIDU '+' (BNI '-' BNP1) ;
DUN = UNSANI '*' RESP1 ;
'SI' VRAI ;
* On the ghost cells we impose DUN = 0
DUN = DUN '*' LIMZER ;
'FINSI' ;
DRN = 'EXCO' 'RN' DUN 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') DUN ('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' DUN 'SCAL' ;
RN = RNL0 '+' DRN ;
GN = GNL0 '+' DGN ;
RETN = RETNL0 '+' DRETN ;
ERRINF = 'MAXIMUM' (DUN '-' DUN0) 'ABS' ;
'SI' (ERRINF < EPSJAC) ;
'MESSAGE'
('CHAINE' 'ITER =' &BL1 ' IJAC =' &BLJAC ' LINF =' ERRINF) ;
'QUITTER' BLJAC ;
'SINON' ;
'SI' (&BLJAC 'EGA' NJAC) ;
'MESSAGE' ('CHAINE' 'ITER =' &BL1
' IJAC =' &BLJAC ' LINF =' ERRINF) ;
'FINSI' ;
'FINSI' ;
'FIN' BLJAC ;
'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
ERRINF = 'MAXIMUM' DUN 'ABS' ;
LISTLINF = LISTLINF 'ET' ('PROG' ERRINF) ;
LISTITER = LISTITER 'ET' ('PROG' &BL1) ;
LISTTPS = LISTTPS 'ET' ('PROG' TPS) ;
'MESSAGE' ;
'MESSAGE'
('CHAINE' 'ITER =' &BL1 ' TPS =' TPS ' LINF =' ERRINF) ;
'MESSAGE' ;
'FINSI' ;
'SI' (((&BL1 '/' FSAF) '*' FSAF) 'EGA' &BL1) ;
'SI' (SAFFAC < SAFFACM) ;
SAFFAC = SAFFAC '*' 2 ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
'FINSI' ;
'FINSI' ;
'SI' (TPS '>EG' TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
'TEMPS' ;
SI GRAPH ;
**********************************************************
**********************************************************
**********************************************************
************** PLOTS *************************************
**********************************************************
**********************************************************
**********************************************************
**********************************************************
*
**** The mesh
*
'TRACER' DOMTOT 'TITR' 'Maillage' ;
*
**** Initial conditions
*
MOD1 = 'MODELISER' ('DOMA' MDOMTOT 'MAILLAGE') 'THERMIQUE' ;
'SI' VRAI ;
VN PN0 = 'PRIM' 'PERFMONO'
RN0 GN0 RETN0 GAMN 'TRICHE' ;
CHM_RN = 'KCHA' MDOMTOT 'CHAM' RN0 ;
CHM_VN = 'KCHA' MDOMTOT 'CHAM' VN ;
CHM_PN = 'KCHA' MDOMTOT 'CHAM' PN0 ;
'TRAC' CHM_RN MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'rho at t=' 0.0) ;
'TRAC' CHM_VN MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'v at t= ' 0.0) ;
'TRAC' CHM_PN MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'p at t= ' 0.0) ;
'FINSI' ;
*
**** The 2D graphics
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN 'TRICHE' ;
'SI' ('EGA' ('TYPE' TPS) 'CHPOINT') ;
tps = 'MINIMUM' TPS ;
'FINSI' ;
'SI' VRAI ;
CHM_RN = 'KCHA' MDOMTOT 'CHAM' RN ;
CHM_VN = 'KCHA' MDOMTOT 'CHAM' VN ;
CHM_PN = 'KCHA' MDOMTOT 'CHAM' PN ;
'TRAC' CHM_RN MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'rho at t=' TPS) ;
'TRAC' CHM_VN MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'v at t= ' TPS) ;
'TRAC' CHM_PN MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'p at t= ' TPS) ;
'FINSI' ;
everr = 'EVOL' 'MANU' 'niter' LISTITER 'Log(Linf)'
(('LOG' LISTLINF) '/' ('LOG' 10.)) ;
'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at ' tps) ;
'FINSI' ;
*
**** Test de convergence
*
AA = 'EXTRAIRE' LISTLINF ('DIME' LISTLINF) ;
'SI' (AA > 1.0D-8) ;
'MESSAGE' 'Probleme en KONV' ;
'ERREUR' 5 ;
'FINSI' ;
*
*** Production d'entropie
*
SN0 = ('LOG' PN0) '-' (GAMN * ('LOG' RN0)) ;
SN0 = 'EXP' SN0 ;
SN = ('LOG' PN) '-' (GAMN * ('LOG' RN)) ;
SN = 'EXP' SN ;
ERR = SN '-' SN0 ;
CH_ERR = 'KCHA' MDOMTOT 'CHAM' ERR ;
'SI' GRAPH ;
'TRAC' CH_ERR MOD1 ('CONTOUR' DOMTOT)
'TITR' ('CHAINE' 'Delta s at t=' TPS) ;
'FINSI' ;
AA = 'MAXIMUM' ERR 'ABS' ;
'SI' (AA > 1.0d-2) ;
'MESSAGE' 'Probleme en KONV' ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer