* fichier : cavitefmm.dgibi ************************************************************************ * Section : Fluides Thermique ************************************************************************ ************************************************************************* * 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' ;