* 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' ; 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 * * **** LIGD * 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 * 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 ; ELP1 ; 'REPETER' BL1 NBL1 ; ELP1 ; ELP1 ; 'FIN' BL1 ; ELP1 ; ELP1 ; 'REPETER' BL1 NBL1 ; ELP1 ; ELP1 ; 'FIN' BL1 ; 'SI' FAUX ; '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' ; RECIN = 0.5 '*' CELL ; REIN = PN '/' (GAMN '-' 1.0) ; RETN = RECIN '+' REIN ; DETR CELL ; DETR RECIN ; DETR REIN ; 'RESPRO' RVN RETN ; 'FINPROC' ; *************************************************************** *************************************************************** *************************************************************** 'NATURE' 'DISCRET' ; 'UY' 0.0 'NATURE' 'DISCRET') ; 'NATURE' 'DISCRET' ; 'NATURE' 'DISCRET' ; GN0 RETN0 = CONS RN0 VN0 PN0 GAMN ; ERRO = 'MAXIMUM' (PN1 '-' PN0) 'ABS' ; 'SI' (ERRO > 1.0D-6) ; 'MESSAGE' 'Problem in the ic file!!!' 'ERREUR' 5 ; 'FINSI' ; * **** Plot of IC * 'SI' FAUX ; 'TRACER' CHM_RN MOD1 'TRACER' CHM_PN MOD1 'TRACER' CHM_VN MOD1 'FINSI' ; MOTRN = 'RN' ; MOTRNX = 'RUX' ; MOTRNY = 'RUY' ; MOTRETN = 'RETN' ; ********************************************************** ********************************************************** ********************************************************** ************** COMPUTATION OF THE SOLUTION *************** ********************************************************** ********************************************************** ********************************************************** ********************************************************** * **** 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 * GRADRN ALRN MCHRNC = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM GRADPN ALPN MCHPNC = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM GRADVN ALVN MCHVNC = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' TYPELIM * Limiters zero into ghost cells 'NATU' 'DISCRET') 'ET' 'NATU' 'DISCRET') 'ET' 'NATU' 'DISCRET') ; * Names of the residuum components TPS = 0.0 ; * **** Temporal loop * 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'Methode = ' METO) ; 'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ; 'MESSAGE' ; 'TEMPS' ; 'REPETER' BL1 NITER ; * **** Primitive variables * 'SI' LOGREC ; GRADRN ALRN = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM GRADPN ALPN = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' TYPELIM GRADVN ALVN = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' TYPELIM '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' ; MDOMTOT RN GRADRN ALRN0 VN GRADVN ALVN0 PN GRADPN ALPN0 GAMN ; 'SINON' ; MDOMTOT RN VN PN GAMN ; 'FINSI' ; MDOMTOT LISTINCO ROF VITF PF GAMF ; * *** Time step * DTPS = DELTAT '*' SAFFAC ; 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' ; 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' ; 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' ; '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 * 'SI' VRAI ; VN PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RETN0 GAMN 'TRICHE' ; 'FINSI' ; * **** The 2D graphics * tps = 'MINIMUM' TPS ; 'FINSI' ; 'SI' VRAI ; 'FINSI' ; (('LOG' LISTLINF) '/' ('LOG' 10.)) ; 'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at ' tps) ; 'FINSI' ; * **** Test de convergence * '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 ; 'SI' GRAPH ; 'FINSI' ; AA = 'MAXIMUM' ERR 'ABS' ; 'SI' (AA > 1.0d-2) ; 'MESSAGE' 'Probleme en KONV' ; 'ERREUR' 5 ; 'FINSI' ; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales