* fichier : sinebum_fmm2.dgibi
*********************************************************************
* *
* SINE-SHAPED BUMP *
* *
* CALCUL DE L ECOULEMENT SUBSONIQUE ISENTROPIQUE STATIONNAIRE DANS *
* UN CANAL *
* *
* Methode implicite sans matrice *
* *
* BECCANTINI A., SFME/LTMF, MAI 2002 *
* *
* JANVIER 2003: *
* Boundary conditions imposed on the border without using ghost *
* cells. *
* Boundary conditions correctly taken into account into the *
* explicit and implicit part. *
* They are not recomputed at each Jacobi iteration. *
* *
*********************************************************************
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 ;
*
**** Creation of MODELS
*
MDOMINT = 'MODELISER' DOMINT 'EULER' ;
MLIGG = 'MODELISER' LIGG 'EULER' ;
MLIGD = 'MODELISER' LIGD 'EULER' ;
*
**** Creation of DOMAINE tables via the MODEL object
*
TDOMINT = 'DOMA' MDOMINT 'VF' ;
TLIGG = 'DOMA' MLIGG 'VF' ;
TLIGD = 'DOMA' MLIGD 'VF' ;
QDOMINT = TDOMINT . 'QUAF' ;
QLIGD = TLIGD . 'QUAF' ;
QLIGG = TLIGG . 'QUAF' ;
'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QLIGG ;
'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QLIGD ;
'SI' GRAPH ;
'TRACER' DOMINT 'TITRE' 'Domaine' ;
'FINSI' ;
***************************************************************
***************************************************************
***************************************************************
************** Initial conditions *****************************
***************************************************************
***************************************************************
***************************************************************
GAMAIR = 1.4 ;
RO = 1.4 ;
P = 1.0 ;
U = 0.5;
***************************************************************
***************************************************************
***************************************************************
************** Boundary Conditions ****************************
***************************************************************
***************************************************************
***************************************************************
* Subsonic inlet.
* We impose the total enthalpy and the entropy.
* The velocity direction is supposed to be 90 degrees
*
HTIN = (GAMAIR '/' (GAMAIR '-' 1.0)) '*' P ;
ECIN = 0.5 '*' RO '*' U '*' U ;
HTIN = HTIN '+' ECIN ;
HTIN = HTIN '/' RO ;
SINL = P '/' (RO '**' GAMAIR) ;
CHPINL = 'MANUEL' 'CHPO' ('DOMA' MLIGG 'CENTRE') 2 'HT' HTIN
'S' sinl ;
* Subsonic outlet
* We impose the pressure
*
CHPOUT = 'MANUEL' 'CHPO' ('DOMA' MLIGD 'CENTRE') 1 'PN' P ;
***************************************************************
***** 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' MDOMINT 'CENTRE') 1 'SCAL' ro
'NATURE' 'DISCRET' ;
VN0 = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2 'UX' u
'UY' 0.0 'NATURE' 'DISCRET' ;
PN0 = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' p
'NATURE' 'DISCRET' ;
GAMN = 'MANU' 'CHPO' ('DOMA' MDOMINT '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' MDOMINT 'MAILLAGE' ) 'THERMIQUE';
'SI' FAUX ;
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN0 ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN0 ;
CHM_VN = 'KCHA' MDOMINT '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' ;
MOTGNX = 'RUX' ;
MOTGNY = 'RUY' ;
MOTGN = 'MOTS' MOTGNX MOTGNY ;
MOTVNX = 'UX' ;
MOTVNY = 'UY' ;
MOTVN = 'MOTS' MOTVNX MOTVNY ;
MOTRETN = 'RETN' ;
MOTPN = 'PN' ;
LISTINCO = 'MOTS' MOTRN MOTGNX MOTGNY MOTRETN ;
LISTPRIM = 'MOTS' MOTRN MOTVNX MOTVNY MOTPN ;
**********************************************************
**********************************************************
**********************************************************
************** 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 = 4000 ;
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 = 20. ;
FSAF = 400 ;
* Jacobi iterations
NJAC = 15 ;
EPSJAC = 1.0D-16 ;
RN = 'COPIER' RN0 ;
GN = 'COPIER' GN0 ;
RETN = 'COPIER' RETN0 ;
*
**** Coeff to compute gradients
* Boundary conditions have to be taken into account
*
SCALBC = 'MANUEL' 'CHPO' (('DOMA' MLIGD 'CENTRE') 'ET'
('DOMA' MLIGG 'CENTRE')) 1 'SCAL' 1.0 ;
VECTBC = 'MANUEL' 'CHPO' (('DOMA' MLIGD 'CENTRE') 'ET'
('DOMA' MLIGG 'CENTRE')) 2 'UX' 1.0 'UY' 1.0 ;
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
GRADRN ALRN MCHRNC = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') RN 'CLIM' SCALBC ;
GRADPN ALPN MCHPNC = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') PN 'CLIM' SCALBC ;
GRADVN ALVN MCHVNC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' TYPELIM
('MOTS' 'UX' 'UY') VN 'CLIM' VECTBC ;
*
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 ;
*
*** Boundary conditions
* Subsonic inlet and subsonic outlet
*
RCHLIM1 RESLIM1 = 'KONV' 'VF' 'PERFMONO' 'CLIM' 'RESI'
MDOMINT MLIGG LISTINCO LISTPRIM
RN VN PN GAMN CHPINL 'INSU' ;
RCHLIM2 RESLIM2 = 'KONV' 'VF' 'PERFMONO' 'CLIM' 'RESI'
MDOMINT MLIGD LISTINCO LISTPRIM
* RN VN PN GAMN CHPOUT 'OUTP' ;
RN VN PN GAMN CHPOUT 'OUTP' ;
* RCHLIM = primitive variables on the boundary
* RCHRESL = contribution of boundary conditions to the residuum
RESLIM = RESLIM1 '+' RESLIM2 ;
RCHLIM = RCHLIM1 '+' RCHLIM2 ;
MAILLIM = 'EXTR' RCHLIM 'MAILLAGE' ;
*
'SI' LOGREC ;
GRADRN ALRN = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') RN
'CLIM' ('EXCO' MOTRN RCHLIM 'SCAL') 'GRADGEO' MCHRNC ;
GRADPN ALPN = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') PN
'CLIM' ('EXCO' MOTPN RCHLIM 'SCAL') 'GRADGEO' MCHPNC ;
GRADVN ALVN = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' TYPELIM
('MOTS' 'UX' 'UY') VN
'CLIM' ('EXCO' MOTVN RCHLIM MOTVN) 'GRADGEO' MCHVNC ;
'SI' (&BL1 < NLCB) ;
* Limiters zero on ghost cells
ALRN0 = 'COPIER' ALRN ;
ALPN0 = 'COPIER' ALPN ;
ALVN0 = 'COPIER' ALVN ;
'SINON' ;
'SI' (&BL1 'EGA' NLCB) ;
'MESSAGE' ;
'MESSAGE' 'On gele les limiteurs!!!' ;
'MESSAGE' ;
'FINSI' ;
'FINSI' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
MDOMINT
RN GRADRN ALRN0
VN GRADVN ALVN0
PN GRADPN ALPN0
GAMN ;
'SINON' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
MDOMINT
RN VN PN GAMN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
MDOMINT LISTINCO ROF VITF PF GAMF MAILLIM ;
RESIDU = RESIDU '+' RESLIM ;
*
*** 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 ;
*
*** We do take into account the boundary conditions into 'AN' and
* 'BN'
ANI = 'KONV' 'VF' 'PMONOFMM' 'AN' LISTINCO
MDOMINT RN GN RETN GAMN SAFFAC 'CLIM' LISTPRIM RCHLIM ;
BNI = 'KONV' 'VF' 'PMONOFMM' 'BN' LISTINCO
MDOMINT RN GN RETN GAMN 'CLIM' LISTPRIM RCHLIM ;
UNSANI = ('INVERSE' ANI) ;
DUN = UNSANI '*' RESIDU ;
DRN = 'EXCO' MOTRN DUN 'SCAL' ;
DGN = 'EXCO' MOTGN DUN MOTVN ;
DRETN = 'EXCO' MOTRETN DUN 'SCAL' ;
RN = RNL0 '+' DRN ;
GN = GNL0 '+' DGN ;
RETN = RETNL0 '+' DRETN ;
'REPETER' BLJAC NJAC ;
DUN0 = DUN ;
BNP1 = 'KONV' 'VF' 'PMONOFMM' 'BN' LISTINCO
MDOMINT RN GN RETN GAMN 'CLIM' LISTPRIM RCHLIM ;
RESP1 = RESIDU '+' (BNI '-' BNP1) ;
DUN = UNSANI '*' RESP1 ;
DRN = 'EXCO' MOTRN DUN 'SCAL' ;
DGN = 'EXCO' MOTGN DUN MOTVN ;
DRETN = 'EXCO' MOTRETN 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' DOMINT 'TITR' 'Maillage' ;
*
**** Initial conditions
*
MOD1 = 'MODELISER' ('DOMA' MDOMINT 'MAILLAGE') 'THERMIQUE' ;
'SI' VRAI ;
VN PN0 = 'PRIM' 'PERFMONO'
RN0 GN0 RETN0 GAMN 'TRICHE' ;
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN0 ;
CHM_VN = 'KCHA' MDOMINT 'CHAM' VN ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN0 ;
'TRAC' CHM_RN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'rho at t=' 0.0) ;
'TRAC' CHM_VN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'v at t= ' 0.0) ;
'TRAC' CHM_PN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'p at t= ' 0.0) ;
'FINSI' ;
*
**** The 2D graphics
*
VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN 'TRICHE' ;
CN2 = GAMN '*' (PN '/' RN) ;
VN2 = 'PSCAL' VN VN MOTVN MOTVN ;
MACHN2 = VN2 '/' CN2 ;
MACHN = MACHN2 '**' 0.5 ;
'SI' ('EGA' ('TYPE' TPS) 'CHPOINT') ;
tps = 'MINIMUM' TPS ;
'FINSI' ;
'SI' VRAI ;
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN ;
CHM_VN = 'KCHA' MDOMINT 'CHAM' VN ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN ;
CHM_MN = 'KCHA' MDOMINT 'CHAM' MACHN ;
'TRAC' CHM_RN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'rho at t=' TPS) ;
'TRAC' CHM_VN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'v at t= ' TPS) ;
'TRAC' CHM_PN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'p at t= ' TPS) ;
'TRAC' CHM_MN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'Mach 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' MDOMINT 'CHAM' ERR ;
'SI' GRAPH ;
'OPTION' 'ISOV' 'SULI' ;
'TRAC' CH_ERR MOD1 ('CONTOUR' DOMINT)
'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