* fichier : cavitefmm.dgibi
*************************************************************************
* CAVITE CARRE A PAROI DEFILANTE *
* Methode implicite sans matrice pour les equations de *
* Navier-Stokes (bas Mach) *
* *
* BECCANTINI A., SFME/LTMF, DEC 2003 *
* *
*************************************************************************
*
* The real file starts after the procedure PNSSM
*
*************************************************************************
*************************************************************************
*********** SOLUTION OF THE NAVIER-STOKES EQUATIONS ********************
*********** FMM ********************
*************************************************************************
*************************************************************************
*
* RVX . 'RESULTS' : table containing the results
*
* RVX . 'FREQI' : frequency of presenting the results
*
* RVX . 'MODEL' : model object
* RVX . 'RN0' : density
* RVX . 'GN0' : qdm
* RVX . 'RET0' : total energy
*
* RVX . 'PGAS' : table containing the gas model
* RVX . 'PGAS' . 'GAMN' : gamma
* RVX . 'PGAS' . 'MU' : dynamic viscosity (kg/m^3 x m^2/s)
* RVX . 'PGAS' . 'LAMBDA' : heat diffusion (W/K/m)
*
*
* RVX . 'GRAVITY' : gravity
*
* RVX . 'LISTCONS' : name of the conservative variables
* RVX . 'LISTPRIM' : name of the primitive variables
*
* RVX . 'LISTERR' : name of the error variables
*
* RVX . 'METHOD' : numerical scheme
* RVX . 'CUTOFF' : cut off speed
*
* RVX . 'SPACEA' : space accuracy
* RVX . 'LIMITER' : limiter type
* RVX . 'TIMEA' : time accuracy
*
* RVX . 'T0' : initial time
* RVX . 'TFINAL' : final time
* RVX . 'DTPS' ('CFL') : time step or CFL number
*
* RVX . 'DCFL' : CFL number for dual time
*
* Error criteria for dual time loop:
* RVX . 'NDTITER' : number of iterations
* RVX . 'RELERR' : Logical (Relative error or absolute error)
* RVX . 'EPSDT' : error at which dual time iterations are stopped
*
* RVX . 'NJAC' : Jacobi iterations
*
* RVX . 'PROLIM' : table called by the procedure PROLIM,
* procedure to compute boundary conditions
*
* RVX . 'DIFTIMP' : boundary condition on temperature
* RVX . 'DIFGTIMP' : boundary condition on gradient of temperature
* RVX . 'DIFVIMP' : boundary condition on speed
* RVX . 'DIFGVIMP' : boundary condition on gradient of speed
* RVX . 'DIFTAUI' : boundary condition on constraint tensor
* RVX . 'DIFQIMP' : boundary condition on heat flux
*
*
'DEBPROC' PNSSM ;
'ARGUMENT' RVX*TABLE ;
*
*
'SI' ('EXISTE' RVX 'RESULTS') ;
'MESS' 'Table RESULTS already exists' ;
'SINON' ;
RVX . 'RESULTS' = 'TABLE' ;
RVX . 'RESULTS' . 'TPS' = RVX . 'T0' ;
RVX . 'RESULTS' . 'RN' = 'COPIER' (RVX . 'RN0') ;
RVX . 'RESULTS' . 'GN' = 'COPIER' (RVX . 'GN0') ;
RVX . 'RESULTS' . 'RET' = 'COPIER' (RVX . 'RET0') ;
RVX . 'RESULTS' . 'LISTLINF' = 'PROG' ;
RVX . 'RESULTS' . 'LISTITDT' = 'LECT' ;
RVX . 'RESULTS' . 'LISTITER' = 'LECT' ;
RVX . 'RESULTS' . 'NITER' = 0 ;
'FINSI' ;
MDOMINT = RVX . 'MODEL' ;
NELT = 'NBEL' ('DOMA' MDOMINT 'CENTRE') ;
*
***** Physical properties
*
GAMN = RVX . 'PGAS' . 'GAMN' ;
GAMSCA1 = 'MAXIMUM' GAMN ;
GAMSCA2 = 'MINIMUM' GAMN ;
'SI' ('EGA' GAMSCA2 GAMSCA1 0.0001) ;
GAMSCAL = GAMSCA1 ;
'SINON' ;
'MESSAGE' ;
'MESSAGE' 'Gamma is not constant' ;
'ERREUR' 21 ;
'FINSI' ;
*
MU = (RVX . 'PGAS' . 'MU') ;
LAMBDA = (RVX . 'PGAS' . 'LAMBDA') ;
R = (RVX . 'PGAS' . 'R') ;
CV = R '/' (gamscal '-' 1.) ;
*
LISTINCO = RVX . 'LISTCONS' ;
LISTPRIM = RVX . 'LISTPRIM' ;
LISTERR = RVX . 'LISTERR' ;
* Names of the gradient of temperature
LMGTEMP = 'MOTS' 'P1DX' 'P1DY' ;
* Names of the gradient of speed
LMGVIT = 'MOTS' 'P1DX' 'P1DY' 'P2DX' 'P2DY' ;
*
* Upwind scheme
*
METO = RVX . 'METHOD' ;
* Space accuracy (1 or 2) and limiter
* Time accuracy (1 or 2)
ORDESP = RVX . 'SPACEA' ;
TYPELIM = RVX . 'LIMITER' ;
ORDTPS = RVX . 'TIMEA' ;
* Initial/final time
* Deltat or CFL
TPS = RVX . 'RESULTS' . 'TPS' ;
TFINAL = RVX . 'TFINAL' ;
'SI' ('EXISTE' RVX 'DTPS') ;
'SI' ('EXISTE' RVX 'CFL') ;
'MESSAGE' 'DTPS or CFL ???' ;
'ERREUR' 21 ;
'FINSI' ;
DTPS = RVX . 'DTPS' ;
'SINON' ;
CFL = RVX . 'CFL' ;
'FINSI' ;
* Dual time iterations
NDT = RVX . 'NDTITER' ;
* Relative error
EPSDT = RVX . 'EPSDT' ;
* Jacobi iterations
* NJAC = RVX . 'NJAC' ;
* Cut off speed
ICO = RVX . 'CUTOFF' ;
* To compute the diffusive cut-off
DELTAX = 'DOMA' MDOMINT 'XXDIEMIN' ;
USDELTAX = 'INVERSE' DELTAX ;
*
**** Conservative variables
*
*
MOT1 = 'EXTRAIRE' LISTINCO 1 ;
'SI' ('EGA' ('VALE' 'DIME') 2) ;
NOMMOM = 'EXTRAIRE' LISTINCO
('LECT' 2 3 ) ;
NOMVEL = 'MOTS' 'UX' 'UY' ;
MOT2 = 'EXTRAIRE' LISTINCO 4 ;
'SINON' ;
NOMMOM = 'EXTRAIRE' LISTINCO
('LECT' 2 3 4) ;
NOMVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
MOT2 = 'EXTRAIRE' LISTINCO 5 ;
'FINSI' ;
RN0 = 'REDU' (RVX . 'RESULTS' . 'RN') ('DOMA' MDOMINT 'CENTRE') ;
GN0 = 'REDU' (RVX . 'RESULTS' . 'GN') ('DOMA' MDOMINT 'CENTRE') ;
RET0 = 'REDU' (RVX . 'RESULTS' . 'RET') ('DOMA' MDOMINT 'CENTRE') ;
*
**** Primitive variables
*
VN0 PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
TN0 = PN0 '/' (R '*' RN0) ;
MOTRN = 'EXTRAIRE' LISTPRIM 1 ;
'SI' ('EGA' ('VALE' 'DIME') 2) ;
MOTVN = 'EXTRAIRE' LISTPRIM
('LECT' 2 3 ) ;
MOTPN = 'EXTRAIRE' LISTPRIM 4 ;
'SINON' ;
MOTVN = 'EXTRAIRE' LISTPRIM
('LECT' 2 3 4) ;
MOTPN = 'EXTRAIRE' LISTPRIM 5 ;
'FINSI' ;
*
**** Coeff to compute gradients for convective term (MCHSCA, MCHVEC)
* Boundary conditions have to be taken into account
*
RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
RN0 VN0 PN0 GAMN ;
MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
'SI' (MAILLIM 'EGA' 0) ;
CHPVID PIPI = 'KOPS' 'MATRIK' ;
MAILLIM = 'DIFF' ('DOMA' MDOMINT 'CENTRE')
('DOMA' MDOMINT 'CENTRE') ;
GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN0 ;
GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
NOMVEL GN0 'CLIM' (RVX . 'DIFVIMP') ;
'SINON' ;
SCALBC = 'MANUEL' 'CHPO' MAILLIM 1 'SCAL' 1.0 ;
'SI' (('VALEUR' 'DIME') 'EGA' 2) ;
VECTBC = 'MANUEL' 'CHPO' MAILLIM
'NATURE' 'DISCRET' 2 'UX' 0.0 'UY' 0.0 ;
* Same name as NOMVEL
'SINON' ;
VECTBC = 'MANUEL' 'CHPO' MAILLIM
'NATURE' 'DISCRET' 3 'UX' 0.0 'UY' 0.0 'UZ' ;
'FINSI' ;
VECTBC = VECTBC 'ET' (RVX . 'DIFVIMP') ;
*
GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN0 'CLIM' SCALBC ;
GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
NOMVEL VN0 'CLIM' VECTBC ;
'FINSI' ;
*
*** Coeff to compute gradient for diffusive flux
*
GRADTN MCHDIT = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
LMGTEMP TN0 (RVX . 'DIFTIMP') (RVX . 'DIFGTIMP') ;
GRADVN MCHDIV = 'PENT' MDOMINT 'FACE' 'DIAMAN2' NOMVEL LMGVIT
VN0 (RVX . 'DIFVIMP') (RVX . 'DIFGVIMP') ;
NOMVEL = 'MOTS' 'UX' 'UY' ;
*
*
*** After each dual time loop, we could display
* the evolution of the error in the dual time loop
* the evolution of the dual time step (the safety factor)
*
LISTLINF = RVX . 'RESULTS' . 'LISTLINF' ;
LISTITDT = RVX . 'RESULTS' . 'LISTITDT' ;
LISTITER = RVX . 'RESULTS' . 'LISTITER' ;
*
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ;
'TEMPS' 'ZERO' ;
************************************************************************
************************************************************************
**** Temporal loop *****************************************************
************************************************************************
************************************************************************
RN_N1M1 = 'COPIER' RN0 ;
GN_N1M1 = 'COPIER' GN0 ;
RET_N1M1 = 'COPIER' RET0 ;
AA = 'DIME' LISTITER ;
'SI' (AA > 0) ;
PTITER = 'EXTRAIRE' LISTITER AA ;
'SINON' ;
PTITER = 0 ;
'FINSI' ;
DUSDT = 0.0D0 ;
'REPETER' BLITER -1 ;
PTITER = PTITER '+' 1 ;
*
**** Personal procedure
*
PROCPT RVX ;
*
*
**** _N1M = (t^n,\tau^m)
* _N1M1 = (t^n,\tau^{m+1})
*
*
************************************************************************
****** Loop on dual time***********************************************
************************************************************************
*
*
*** DUSDT0 is the increment of DUSDT in the previous (physical) time
* iteration.
*
DUSDT0 = DUSDT ;
DUSDT = 0.0D0 ;
'REPETER' BLDT NDT ;
RN_N1M = RN_N1M1 ;
GN_N1M = GN_N1M1 ;
RET_N1M = RET_N1M1 ;
*
**** Primitive variables
*
VN_N1M PN_N1M = 'PRIM' 'PERFMONO' RN_N1M GN_N1M RET_N1M GAMN ;
*
*** Boundary conditions
*
RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
RN_N1M VN_N1M PN_N1M GAMN ;
*
****** First/second order reconstruction
*
'SI' (ORDESP 'EGA' 2) ;
*
NNLIM = 'NBNO' MAILLIM ;
'SI' (NNLIM 'EGA' 0) ;
RNLIM = CHPVID ;
PNLIM = CHPVID ;
VNLIM = RVX . 'DIFVIMP' ;
'SINON' ;
RNLIM = ('EXCO' MOTRN RCHLIM 'SCAL') ;
PNLIM = ('EXCO' MOTPN RCHLIM 'SCAL') ;
VNLIM = ('EXCO' MOTVN RCHLIM MOTVN)
'ET' (RVX . 'DIFVIMP') ;
'FINSI' ;
GRADRN ALRN0 = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') RN_N1M
'CLIM' RNLIM 'GRADGEO' MCHSCA ;
GRADPN ALPN0 = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') PN_N1M
'CLIM' PNLIM 'GRADGEO' MCHSCA ;
GRADVN ALVN0 = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' TYPELIM
('MOTS' 'UX' 'UY') VN_N1M
'CLIM' VNLIM 'GRADGEO' MCHVEC ;
* 'SI' (&BLDT < NLCB) ;
* ALRN0 = 'COPIER' ALRN ;
* ALPN0 = 'COPIER' ALPN ;
* ALVN0 = 'COPIER' ALVN ;
* 'SINON' ;
* 'SI' (&BLDT 'EGA' NLCB) ;
* 'MESSAGE' ;
* 'MESSAGE' 'On gele les limiteurs!!!' ;
* 'MESSAGE' ;
* 'FINSI' ;
* 'FINSI' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
MDOMINT
RN_N1M GRADRN ALRN0
VN_N1M GRADVN ALVN0
PN_N1M GRADPN ALPN0
GAMN ;
'SINON' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
MDOMINT
RN_N1M VN_N1M PN_N1M GAMN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
MDOMINT LISTINCO ROF VITF PF GAMF MAILLIM
ICO (MU '*' ('INVERSE' RN_N1M) '*' USDELTAX) ;
* ICO ICO ;
RESIDU = RESIDU '+' RESLIM ;
*
**** La gravite
*
RESGRA = 'FIMP' 'VF' 'GRAVMONO' 'RESI' LISTINCO
RN_N1M GN_N1M1 (RVX . 'GRAVITY') ;
RESIDU = RESIDU '+' RESGRA ;
*
**************************************
**** Diffusive terms *****************
**************************************
*
TN_N1M = PN_N1M '/' (R '*' RN_N1M) ;
*
GRADTN = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
LMGTEMP TN_N1M (RVX . 'DIFTIMP') (RVX . 'DIFGTIMP')
'GRADGEO' MCHDIT ;
GRADVN = 'PENT' MDOMINT 'FACE' 'DIAMAN2' NOMVEL LMGVIT
VN_N1M (RVX . 'DIFVIMP') (RVX . 'DIFGVIMP')
'GRADGEO' MCHDIV ;
NOMVEL = 'MOTS' 'UX' 'UY' ;
*
ICACCA RESIDI DTCACCA = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'EXPL'
MDOMINT MU LAMBDA CV RN_N1M VN_N1M TN_N1M GRADVN GRADTN
LISTINCO 'VIMP' (RVX . 'DIFVIMP') 'TAUI' (RVX . 'DIFTAUI')
'QIMP' (RVX . 'DIFQIMP') ;
RESIDU = RESIDU '+' RESIDI ;
COEFV = (0.0 '*' RN_N1M) '+' 1.0 ;
COEFV = COEFV '*' (lambda '*' (gamscal '-' 1.)) ;
COEFV = COEFV '/' (R '*' RN_N1M) ;
* 'MESSAGE' ;
* 'MESSAGE' 'Je fais n importe' ;
* COEFV = COEFV '*' 0.25 ;
*
****** Residuum for dual tims stepping also involved the
* variation of the conserved variables with respect
* to time
'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
RESIDU = RESIDU '-' DUSDT ;
'SINON' ;
RESIDU = RESIDU '-' ((1.5 '*' DUSDT) '-' (0.5 '*' DUSDT0)) ;
'FINSI' ;
*
*** Time step at the first iteration/jacobi iteration
*
'SI' (&BLDT 'EGA' 1) ;
'SI' ('EXISTE' RVX 'CFL') ;
DTPS = (RVX . 'CFL') '*' 2.0D0 '*' DELTAT ;
'SINON' ;
DTPS = RVX . 'DTPS' ;
'FINSI' ;
DTPS = 'MINIMUM' ('PROG' DTPS ((TFINAL '-' TPS) '*' 1.001)) ;
TPS = TPS '+' DTPS ;
*
NJAC = 'ENTIER' ('MINIMUM' (RVX . 'NJACITER')) ;
NJAC0 = NJAC ;
'SINON' ;
NJAC = 'IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
(RVX . 'NJACLERR') (RVX . 'NJACITER') ;
NJAC = 'MAXIMUM' ('PROG' NJAC0 NJAC) ;
NJAC = 'ENTIER' NJAC ;
NJAC0 = NJAC ;
'FINSI' ;
*
*** JACOBI
*
*
**** CFL dual
*
'SI' (&BLDT 'EGA' 1) ;
SAFFACD = ('MINIMUM' (RVX . 'DCFL')) '*' 2 ;
'SINON' ;
SAFFACD = ('IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
(RVX . 'DCFLERR') (RVX . 'DCFL')) '*' 2 ;
'FINSI' ;
'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
* DUN IPRO = 'DETO' (RVX . 'TYPEJAC')
DUN IPRO = 'KONV' 'VF' 'PMON1FMM' (RVX . 'TYPEJAC')
LISTINCO MDOMINT RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
DTPS SAFFACD
NJAC 'CLIM' LISTPRIM RCHLIM COEFV ;
'SINON' ;
* DUN IPRO = 'DETO' (RVX . 'TYPEJAC')
DUN IPRO = 'KONV' 'VF' 'PMON1FMM' (RVX . 'TYPEJAC')
LISTINCO MDOMINT RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
(DTPS '/' 1.5)
SAFFACD NJAC 'CLIM' LISTPRIM RCHLIM COEFV ;
'FINSI' ;
* 'FINSI' ;
'SI' (IPRO 'NEG' 0) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dans FMM' ;
'MESSAGE' ;
'ERREUR' 21 ;
'FINSI' ;
*
**** We compute DUSDT for the future loop
*
DUSDT = DUSDT '+' (DUN '/' DTPS) ;
*
**** We evaluate the conservative variables at t^{n+1}, \tau^{m+1}
*
DRN = 'EXCO' MOT1 DUN 'SCAL' ;
DGN = 'EXCO' NOMMOM DUN NOMVEL ;
DRET = 'EXCO' MOT2 DUN 'SCAL' ;
RN_N1M1 = RN_N1M '+' DRN ;
GN_N1M1 = GN_N1M '+' DGN ;
RET_N1M1 = RET_N1M '+' DRET ;
ERRINF = 'MAXIMUM' DUN 'ABS' LISTERR ;
LISTLINF = LISTLINF 'ET' ('PROG' ERRINF) ;
LISTITDT = LISTITDT 'ET' ('LECT' &BLDT) ;
LISTITER = LISTITER 'ET' ('LECT' PTITER) ;
'SI' ((&BLDT 'EGA' 1) 'OU'
(((&BLDT '/' (RVX . 'FREQI')) '*' (RVX . 'FREQI'))
'EGA' &BLDT)) ;
'MESSAGE' ;
'MESSAGE'
('CHAINE' 'ITER =' PTITER ' TPS =' TPS
' DTITER =' &BLDT ' LINF =' ERRINF
' DCFL =' SAFFACD ' NJAC =' NJAC) ;
'MESSAGE' ;
'FINSI' ;
*
*
*** Update of RVX . 'RESULTS'
*
RVX . 'RESULTS' . 'RN' = RN_N1M1 ;
RVX . 'RESULTS' . 'GN' = GN_N1M1 ;
RVX . 'RESULTS' . 'RET' = RET_N1M1 ;
RVX . 'RESULTS' . 'LISTLINF' = LISTLINF ;
RVX . 'RESULTS' . 'LISTITDT' = LISTITDT ;
RVX . 'RESULTS' . 'LISTITER' = LISTITER ;
*
'SI' (RVX . 'RELERR') ;
* Relative error
'SI' (&BLDT 'EGA' 1) ;
ERRINF0 = ERRINF ;
'SINON' ;
'SI' (ERRINF < (EPSDT '*' ERRINF0)) ;
'QUITTER' BLDT ;
'FINSI' ;
'FINSI' ;
'SINON' ;
'SI' (ERRINF < EPSDT) ;
'QUITTER' BLDT ;
'FINSI' ;
'FINSI' ;
'FIN' BLDT ;
************************************************************************
****** End of the loop on dual time*************************************
************************************************************************
*
*** Update of RVX . 'RESULTS'
*
RVX . 'RESULTS' . 'TPS' = TPS ;
RVX . 'RESULTS' . 'NITER' = (RVX . 'RESULTS' . 'NITER') '+' 1 ;
*
'SI' (TPS '>EG' TFINAL) ;
'QUITTER' BLITER ;
'FINSI' ;
'FIN' BLITER ;
TCPU = 'TEMPS' 'NOEC' ;
RVX . 'RESULTS' . 'TCPU' = TCPU ;
'FINPROC' ;
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
'OPTION' 'DIME' 2
'ELEM' 'QUA4'
'ISOV' 'LIGN'
'TRAC' 'PSC'
'ECHO' 1
;
LOGTRI = VRAI ;
* GRAPH = VRAI ;
GRAPH = FAUX ;
*******************************
**** MESH *****************
*******************************
RAF = 8 ;
L = 1 ;
H = L ;
DENINI = (L '/' (8. '*' RAF)) ;
DENCEN = 8 '*' DENINI ;
A1 = 0.0 0.0 ;
A12 = (L '/' 2.) 0.0 ;
A2 = L 0.0 ;
A23 = L (H '/' 2.0) ;
A3 = L H ;
A34 = (L '/' 2.0) H ;
A4 = 0.0 H ;
A41 = 0.0 (H '/' 2.0) ;
A1A2 = (A1 'DROIT' A12 'DINI' DENINI 'DFIN' DENCEN) 'ET'
(A12 'DROIT' A2 'DINI' DENCEN 'DFIN' DENINI);
A2A3 = (A2 'DROIT' A23 'DINI' DENINI 'DFIN' DENCEN) 'ET'
(A23 'DROIT' A3 'DINI' DENCEN 'DFIN' DENINI);
A3A4 = (A3 'DROIT' A34 'DINI' DENINI 'DFIN' DENCEN) 'ET'
(A34 'DROIT' A4 'DINI' DENCEN 'DFIN' DENINI);
A4A1 = (A4 'DROIT' A41 'DINI' DENINI 'DFIN' DENCEN) 'ET'
(A41 'DROIT' A1 'DINI' DENCEN 'DFIN' DENINI);
DOMINT = 'DALLER' A1A2 A2A3 A3A4 A4A1 'PLAN' ;
'SI' LOGTRI ;
'OPTION' 'ELEM' 'TRI3' ;
DOMINT = 'SURFACE' ( A1A2 'ET' A2A3 'ET' A3A4 'ET' A4A1 ) 'PLAN' ;
'FINSI' ;
MDOMINT = 'MODELISER' DOMINT 'EULER' ;
MA1A2 = 'MODELISER' A1A2 'EULER' ;
MA2A3 = 'MODELISER' A2A3 'EULER' ;
MA3A4 = 'MODELISER' A3A4 'EULER' ;
MA4A1 = 'MODELISER' A4A1 'EULER' ;
$DOMINT = 'DOMA' MDOMINT 'VF' ;
$A1A2 = 'DOMA' MA1A2 'VF' ;
$A2A3 = 'DOMA' MA2A3 'VF' ;
$A3A4 = 'DOMA' MA3A4 'VF' ;
$A4A1 = 'DOMA' MA4A1 'VF' ;
QDOMINT = 'DOMA' MDOMINT 'QUAF' ;
QA1A2 = 'DOMA' MA1A2 'QUAF' ;
QA2A3 = 'DOMA' MA2A3 'QUAF' ;
QA3A4 = 'DOMA' MA3A4 'QUAF' ;
QA4A1 = 'DOMA' MA4A1 'QUAF' ;
'ELIMINATION' QDOMINT (DENINI '/' 100) QA1A2 ;
'ELIMINATION' QDOMINT (DENINI '/' 100) QA2A3 ;
'ELIMINATION' QDOMINT (DENINI '/' 100) QA3A4 ;
'ELIMINATION' QDOMINT (DENINI '/' 100) QA4A1 ;
'SI' GRAPH ;
'TRACER' DOMINT 'TITRE' ('CHAINE' 'NBEL =' ('NBEL' DOMINT)) ;
'FINSI' ;
***************************************
******** PHYSICAL PARAMETERS *********
***************************************
gamscal = 1.4 ;
RAIR = 288. ;
*
ro_0 = 1.0 ;
p_0 = 1.0D5 ;
ret_0 = (1. '/' (gamscal '-' 1.)) '*' p_0 ;
CSON = (gamscal * p_0 '/' ro_0) '**' 0.5 ;
*
VParois = 0.01 '*' cson ;
*
Rey = 1000 ;
mu = (ro_0 '*' L '*' VParois) '/' Rey ;
*
Pr = 0.72 ;
* lambda = mu * cp '/' pr ;
lambda = mu '*' (gamscal '*' RAIR '/' (gamscal '-' 1.0))
'/' Pr ;
***************************************
**** INITIAL CONDITIONS ***************
***************************************
*
* A4 -> A3
*
*
*
*
* A1 A2
*
GAMN = 'KCHT' $DOMINT 'SCAL' 'CENTRE' gamscal ;
RN0 = 'KCHT' $DOMINT 'SCAL' 'CENTRE' ro_0 ;
GN0 = 'KCHT' $DOMINT 'VECT' 'CENTRE' ((ro_0 '*' 0) 0.0) ;
RETN0 = 'KCHT' $DOMINT 'SCAL' 'CENTRE' (ret_0) ;
VN0 PN0 = 'PRIM' 'PERFMONO'
RN0 GN0 RETN0 GAMN ;
TN0 = PN0 '/' (Rair '*' RN0) ;
CHPVID CACCA = 'KOPS' MATRIK ;
* Conditions aux bords
VLIM = ('MANUEL' 'CHPO' ($A3A4 . 'CENTRE') 2 'UX' Vparois 'UY' 0.0
'NATU' 'DISCRET' ) 'ET'
('MANUEL' 'CHPO' (($A4A1 . 'CENTRE') 'ET'
($A1A2 . 'CENTRE') 'ET' ($A2A3 . 'CENTRE'))
2 'UX' 0.0 'UY' 0.0 'NATU' 'DISCRET') ;
GRADVLIM = CHPVID ;
TAULIM = 'COPIER' CHPVID ;
*
TLIM = CHPVID ;
GRADTLIM = 'MANUEL' 'CHPO' (($A3A4 . 'CENTRE') 'ET' ($A4A1 . 'CENTRE')
'ET' ($A1A2 . 'CENTRE') 'ET' ($A2A3 . 'CENTRE'))
2 'P1DX' 0.0 'P1DY' 0.0 'NATU' 'DISCRET' ;
QLIM = 'MANUEL' 'CHPO' (($A3A4 . 'CENTRE') 'ET' ($A4A1 . 'CENTRE')
'ET' ($A1A2 . 'CENTRE') 'ET' ($A2A3 . 'CENTRE'))
2 'UX' 0.0 'UY' 0.0 'NATU' 'DISCRET' ;
*
***** Pictures of the initial conditions
*
MOD1 = 'MODELISER' ($DOMINT . 'MAILLAGE') 'THERMIQUE' ;
'SI' GRAPH ;
CELL = 'VALEUR' 'ISOV' ;
'OPTION' 'ISOV' 'SULI' ;
*
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN0 ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN0 ;
CHM_TN = 'KCHA' MDOMINT 'CHAM' TN0 ;
CHM_VN = 'KCHA' MDOMINT 'CHAM' VN0 ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' 0.0);
TRAC CHM_PN MOD1 'TITR' ('CHAINE' 'P at t=' 0.0);
TRAC CHM_TN MOD1 'TITR' ('CHAINE' 'T at t=' 0.0);
TRAC CHM_VN MOD1 'TITR' ('CHAINE' 'VN at t=' 0.0);
*
'OPTION' 'ISOV' CELL ;
'FINSI' ;
* Names of conserved variables
NOMDEN = 'MOTS' 'RN' ;
NOMMOM = 'MOTS' 'RUX' 'RUY' ;
NOMRET = 'MOTS' 'RETN' ;
NOMVEL = 'MOTS' 'UX' 'UY' ;
NOMPRE = 'MOTS' 'NOMPRE' ;
LISTCONS = NOMDEN 'ET' NOMMOM 'ET' NOMRET ;
LISTP = NOMDEN 'ET' NOMVEL 'ET' NOMPRE ;
UNCONS = ('NOMC' ('MOTS' 'SCAL') NOMDEN RN0 'NATU' 'DISCRET') 'ET'
('NOMC' NOMVEL GN0 NOMMOM 'NATU' 'DISCRET') 'ET'
('NOMC' ('MOTS' 'SCAL') RETN0 NOMRET 'NATU' 'DISCRET') ;
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
************ COMPUTATION OF THE SOLUTION *****************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
******PROCEDURE PROCPT **************************************************
*************************************************************************
*
*
'DEBPROC' PROCPT ;
'FINPROC' ;
*
*************************************************************************
******PROCEDURE PROLIM **************************************************
*************************************************************************
*
'DEBPROC' PROLIM ;
'ARGUMENT' RVX*'TABLE' MDOMINT*'MMODEL' LISTINCO*'LISTMOTS'
LISTPRIM*'LISTMOTS'
RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
*
*** Boundary conditions
*
CHPVID PIPI = 'KOPS' MATRIK ;
* 'RESPRO' RCHLIM RESLIM ;
'RESPRO' CHPVID CHPVID ;
'FINPROC' ;
*************************************************************************
****** FIN PROCEDURE PROLIM *********************************************
*************************************************************************
*
RV = 'TABLE' ;
RV . 'FREQI' = 1 ;
RV . 'MODEL' = MDOMINT ;
*
**** Conservative variables / primitive variables
*
*
RV . 'LISTCONS' = LISTCONS ;
RV . 'LISTPRIM' = LISTP ;
*
RV . 'RN0' = RN0 ;
RV . 'GN0' = GN0 ;
RV . 'RET0' = RETN0 ;
*
**** Gas property/gravity
*
RV . 'PGAS' = 'TABLE' ;
RV . 'PGAS' . 'GAMN' = GAMN ;
RV . 'PGAS' . 'MU' = mu ;
RV . 'PGAS' . 'R' = Rair ;
RV . 'PGAS' . 'LAMBDA' = lambda ;
*
RV . 'GRAVITY' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2 'UX' 0.0
'UY' 0.0 ;
*
* Table for BC (convective flux)
*
RV . 'PROLIM' = 'TABLE' ;
*
* BC (diffusive flux)
*
* RVX . 'DIFTIMP' : boundary condition on temperature
* RVX . 'DIFGTIMP' : boundary condition on gradient of temperature
* RVX . 'DIFVIMP' : boundary condition on speed
* RVX . 'DIFGVIMP' : boundary condition on gradient of speed
* RVX . 'DIFTAUI' : boundary condition on constraint tensor
* RVX . 'DIFQIMP' : boundary condition on heat flux
*
*
RV . 'DIFTIMP' = TLIM ;
RV . 'DIFGTIMP' = GRADTLIM ;
RV . 'DIFQIMP' = QLIM ;
RV . 'DIFVIMP' = VLIM ;
RV . 'DIFGVIMP' = GRADVLIM ;
RV . 'DIFTAUI' = TAULIM ;
* GRADTF MCHTNF = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
* ('MOTS' 'P1DX' 'P1DY') TN0 TLIM GRADTLIM ;
* GRADVF MCHVNF = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'UX' 'UY')
* ('MOTS' 'P1DX' 'P1DY' 'P2DX' 'P2DY') VN0 VLIM TAULIM ;
*
* ICACCA RESIDI DTCACCA = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'EXPL'
* MDOMINT MU LAMBDA CV RN0 VN0 TN0 GRADVF GRADTF
* (RV . 'LISTCONS') 'VIMP' (RV . 'DIFVIMP') 'TAUI' (RV . 'DIFTAUI')
* 'QIMP' (RV . 'DIFQIMP') ;
*
**** Numerical parameters
*
*
* Variable to compute Linf error
*
RV . 'LISTERR' = NOMMOM ;
* RV . 'LISTERR' = NOMDEN ;
* RV . 'LISTERR' = NOMRET ;
*
* Upwind scheme
*
RV . 'METHOD' = 'RUSANOLM' ;
* RV . 'METHOD' = 'AUSMPLM' ;
*
**** Cut off
* We take a constant cutoff
*
RV . 'CUTOFF' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1
'SCAL' Vparois ;
* Reconstruction/limiter
* Time accuracy (1 or 2)
* Iterations
* Final time
RV . 'SPACEA' = 2 ;
RV . 'LIMITER' = 'NOLIMITE' ;
* RV . 'LIMITER' = 'LIMITEUR' ;
RV . 'TIMEA' = 1 ;
*
**** Phisical time
*
RV . 'T0' = 0 ;
RV . 'TFINAL' = 1.0D6 ;
RV . 'CFL' = 1.0D16 ;
* RV . 'DTPS' = 1.0D16 ;
*
**** Dual time
*
* Safety factor for the dual time step
* Max. Dual time iterations
* Absolute/relative error
*
RV . 'DCFLERR' = 'PROG' 16. -32. ;
RV . 'DCFL' = 'PROG' 1.0D1 1.0D1 ;
RV . 'NDTITER' = 1000 ;
RV . 'RELERR' = FAUX ;
RV . 'EPSDT' = 1.0D-5 ;
**** Jacobi iterations
RV . 'TYPEJAC' = 'PJACO' ;
* RV . 'TYPEJAC' = 'LJACOF' ;
* RV . 'TYPEJAC' = 'LJACOB' ;
* RV . 'TYPEJAC' = 'LJACOFB' ;
RV . 'NJACITER' = 'PROG' 15 15 ;
RV . 'NJACLERR' = 'PROG' 16. -16 ;
*
**** Parameters for PROCPT
*
'TEMPS' 'ZERO' ;
PNSSM RV ;
* 'TEMPS' 'IMPR' ;
**********************************************************
**********************************************************
**********************************************************
************** PLOTS *************************************
**********************************************************
**********************************************************
**********************************************************
**********************************************************
RN = RV . 'RESULTS' . 'RN' ;
GN = RV . 'RESULTS' . 'GN' ;
RET = RV . 'RESULTS' . 'RET' ;
LISTITER = RV . 'RESULTS' . 'LISTITER' ;
LISTITDT = RV . 'RESULTS' . 'LISTITDT' ;
LISTLINF = RV . 'RESULTS' . 'LISTLINF' ;
*
*** Convergence evolution inside of each iteration
*
NITERE = 'DIME' LISTITER ;
I1 = 1 ;
AA = 'PROG' ;
BB = 'PROG' ;
CC = 'PROG' ;
'REPETER' BLITER NITERE ;
I0 = I1 ;
I1 = 'EXTRAIRE' LISTITER &BLITER ;
'SI' (I1 'EGA' I0) ;
AA = AA 'ET' ('PROG' ('EXTRAIRE' LISTITDT &BLITER)) ;
BB = BB 'ET' ('PROG' ('EXTRAIRE' LISTLINF &BLITER)) ;
'SINON' ;
everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
(('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
('LOG' 10.)) ;
'SI' GRAPH ;
'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
(I1 '-' 1)) ;
'FINSI' ;
AA = 'PROG' ;
BB = 'PROG' ;
'FINSI' ;
'FIN' BLITER ;
everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
(('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
('LOG' 10.)) ;
'SI' GRAPH ;
'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
(I1)) ;
'FINSI' ;
'FINSI' ;
VN PN = 'PRIM' 'PERFMONO' RN GN RET GAMN 'TRICHE' ;
CN2 = GAMN '*' (PN '/' RN) ;
VN2 = 'PSCAL' VN VN NOMVEL NOMVEL ;
MACHN2 = VN2 '/' CN2 ;
MACHN = MACHN2 '**' 0.5 ;
HTN = (GAMN '/' (GAMN '-' 1.0)) '*' (PN '/' RN) ;
ECIN = 0.5 '*' ('PSCAL' VN VN NOMVEL NOMVEL) ;
HTN = HTN '+' ECIN ;
SN = PN '/' (RN '**' GAMSCAL) ;
tps = RV . 'RESULTS' . 'TPS' ;
'SI' GRAPH ;
'OPTION' 'ISOV' 'SURF' ;
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN ;
CHM_VN = 'KCHA' MDOMINT 'CHAM' VN ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN ;
CHM_MN = 'KCHA' MDOMINT 'CHAM' MACHN ;
CHM_HTN = 'KCHA' MDOMINT 'CHAM' HTN ;
CHM_SN = 'KCHA' MDOMINT 'CHAM' SN ;
'TRAC' CHM_RN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'rho at t=' TPS) ;
'TRAC' CHM_VN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'v at t= ' TPS) ;
'TRAC' CHM_PN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'p at t= ' TPS) ;
'TRAC' CHM_MN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'Mach at t= ' TPS) ;
'TRAC' CHM_HTN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'ht at t= ' TPS ' hl =' HTL) ;
'TRAC' CHM_SN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'sn at t= ' TPS ' hl =' SL) ;
* 'OPTION' 'ISOV' 'LIGN' ;
* RNV = 'ELNO' $DOMINT ('KCHT' $DOMINT 'SCAL' 'CENTRE' RN) ;
* 'TRACER' DOMINT RNV ('CONTOUR' DOMINT) 15 'TITRE' 'ro';
* PNV = 'ELNO' $DOMINT ('KCHT' $DOMINT 'SCAL' 'CENTRE' PN) ;
* 'TRACER' DOMINT PNV ('CONTOUR' DOMINT) 15 'TITRE' 'p' ;
* MNV = 'ELNO' $DOMINT ('KCHT' $DOMINT 'SCAL' 'CENTRE' MACHN) ;
* 'TRACER' DOMINT MNV ('CONTOUR' DOMINT) 15 'TITRE' 'Mach' ;
* RUX = 'EXCO' GN 'UX' ;
* RUXV = 'ELNO' $DOMINT ('REDU' RUX $DOMINT . 'CENTRE') ;
* 'TRACER' DOMINT RUXV ('CONTOUR' DOMINT) 15
* 'TITRE' 'rux' ;
VECN = 'VECTEUR' VN ;
'TRACER' VECN DOMINT ;
MODVN = 'PSCAL' VN VN NOMVEL NOMVEL ;
MODVN = MODVN '**' 0.5 ;
MODVN = MODVN '+' (('MAXIMUM' MODVN) '*' 0.000001) ;
VN1 = VN '/' MODVN ;
VECN1 = 'VECTEUR' 0.01 VN1 'JAUNE' ;
'TRACER' DOMINT VECN1 ('CONTOUR' DOMINT)
'TITRE' ('CHAINE' 'Normalized VN : ' TITOLO) ;
'FINSI' ;
*
**** Evolution objects
*
TAB1 = 'TABLE' ;
TAB1 . 'TITRE'= 'TABLE' ;
TAB1 . 1 = 'MARQ CROI NOLI';
TAB1 . 'TITRE' . 1 = 'Reference (Su)' ;
TAB1 . 2 = 'REGU' ;
TAB1 . 'TITRE' . 2 = 'Numerical res.' ;
LXSU = 'PROG' 0 0.0625 0.0703 0.0781 0.0938 0.1563 0.2266
0.2344 0.5 0.8047 0.8594 0.9063 0.9453 0.9531 0.9609
0.9688 1 ;
LUY_X_SU = 'PROG' 0 0.2763 0.2918 0.3053 0.3282 0.3711 0.3279
0.3193 0.0243 -0.317 -0.4245 -0.5182 -0.3972 -0.3421 -0.2816
-0.2175 0 ;
EV_UYX = 'EVOL' 'MANU' 'x' LXSU 'uy' LUY_X_SU ;
LYSU = 'PROG' 0. 0.0547 0.0625 0.0703 0.1016 0.1719 0.2812 0.4531
0.5 0.6172 0.7344 0.8516 0.9531 0.9609 0.9688 0.9766 1 ;
LUX_Y_SU = 'PROG' 0 -0.1788 -0.1999 -0.2204 -0.2972 -0.3834 -0.2759
-0.1058 -0.0605 0.0564 0.1857 0.3316 0.466 0.5109 0.5743 0.6582 1 ;
EV_UXY = 'EVOL' 'MANU' 'ux' LUX_Y_SU 'y' LYSU ;
'OPTION' 'ELEM' 'QUA4' ;
LISTX = 'PROG' 0.5 ;
'REPETER' BL1 ('DIME' LISTX) ;
XL = 'EXTRAIRE' LISTX &BL1 ;
YMIN = 'MINIMUM' ('COORDONNEE' 2 DOMINT) ;
YMAX = 'MAXIMUM' ('COORDONNEE' 2 DOMINT) ;
DELTAY = YMAX '-' YMIN ;
PL1 = XL (YMIN '+' (1.0D-10 '*' DELTAY)) ;
PL2 = XL 0.0 ;
PL3 = XL (YMAX '-' (1.0D-10 '*' DELTAY)) ;
LIG1 = (PL1 'DROIT' PL2 'DINI' DENINI 'DFIN' DENCEN)
'ET' (PL2 'DROIT' PL3 'DINI' DENCEN 'DFIN' DENINI) ;
'SI' GRAPH ;
'TRACER' (DOMINT 'ET' (LIG1 'COULEUR' 'ROUG')) 'TITR' 'lig1' ;
'FINSI' ;
UX = ('EXCO' 'UX' VN) '/' Vparois;
UXV = 'ELNO' UX $DOMINT ;
CH_UX = 'CHANGER' 'CHAM' UXV ('DOMA' MDOMINT 'MAILLAGE') ;
UXLIG1 = 'PROI' CH_UX LIG1 ;
EVUX = 'EVOL' 'CHPO' UXLIG1 LIG1 ;
EVUX1 = 'EVOL' 'MANU' 'ux' ('EXTRAIRE' EVUX 'ORDO') 'y'
('EXTRAIRE' EVUX 'ABSC') ;
'SI' GRAPH ;
'DESSIN' (EV_UXY 'ET' EVUX1) 'TITRE'
('CHAINE' 'Nondimensional ux at x = ' XL) 'LEGE'
TAB1 ;
'FINSI' ;
'FIN' BL1 ;
'OPTION' 'ELEM' 'QUA4' ;
LISTY = 'PROG' 0.5 ;
'REPETER' BL1 ('DIME' LISTY) ;
YL = 'EXTRAIRE' LISTY &BL1 ;
XMIN = 'MINIMUM' ('COORDONNEE' 2 DOMINT) ;
XMAX = 'MAXIMUM' ('COORDONNEE' 2 DOMINT) ;
DELTAX = XMAX '-' XMIN ;
PL1 = (XMIN '+' (1.0D-10 '*' DELTAX)) YL;
PL2 = 0.0 YL ;
PL3 = (XMAX '-' (1.0D-10 '*' DELTAX)) YL ;
LIG1 = (PL1 'DROIT' PL2 'DINI' DENINI 'DFIN' DENCEN)
'ET' (PL2 'DROIT' PL3 'DINI' DENCEN 'DFIN' DENINI) ;
'SI' GRAPH ;
'TRACER' (DOMINT 'ET' (LIG1 'COULEUR' 'ROUG')) 'TITR' 'lig1' ;
'FINSI' ;
UY = ('EXCO' 'UY' VN) '/' Vparois;
UYV = 'ELNO' UY $DOMINT ;
CH_UY = 'CHANGER' 'CHAM' UYV ('DOMA' MDOMINT 'MAILLAGE') ;
UYLIG1 = 'PROI' CH_UY LIG1 ;
EVUY = 'EVOL' 'CHPO' UYLIG1 LIG1 ;
'SI' GRAPH ;
'DESSIN' (EV_UYX 'ET' EVUY)
'TITRE' ('CHAINE' 'Nondimensional uy at y = ' YL)
'LEGE' TAB1 ;
'FINSI' ;
'FIN' BL1 ;
* TEST
LERR = 'EXTRAIRE' EVERR 'ORDO' ;
AA = 10. '**' ('EXTRAIRE' LERR 1) ;
BB = 10. '**' ('EXTRAIRE' LERR ('DIME' LERR)) ;
'SI' ((BB '/' AA) '>' 1.0D-5) ;
'ERREUR' 'Probleme de convergence' ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer