* fichier : carre_therper_scalpass.dgibi ************************************************************************ * Section : Thermique Conduction ************************************************************************ ****************************************************************** * CALCUL DU QUARRE A CHOC * * Modele: gaz mono-espece "thermally perfect" * * Gaz mono-espece "thermally perfect" * * Splitting de scalaire passifs * * Controle de la symetrie * * * * A. BECCANTINI, DRN/DMT/SEMT/LTMF, FEVRIER 2000 * *Modif, 10/07/01, syntaxe de PENT changée * *Modif, 25/09/02, syntaxe de KONV changée * ****************************************************************** 'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'ISOV' 'SULI' 'ECHO' 1 'TRAC' 'X'; GRAPH = FAUX ; * GRAPH = VRAI ; * *** Methodes possibles : * * 'VLH', 'SS' * METO = 'VLH' ; * SAFFAC = safety facton on the CFL SAFFAC = 1.0 ; * LOGSO = second order reconstruction LOGSO = VRAI ; TFINAL = 10. ; NITER = 40 ; ZERO = 1.0D-12 ; ************ * MAILLAGE * ************ PCEN = 0.5 0.5 ; A1 = 0.75 0.5 ; A2 = 1.0 0.5 ; A3 = 1.0 0.75 ; A4 = 1.0 1.0 ; A5 = 0.75 1.0 ; A6 = 0.5 1.0 ; A7 = 0.5 0.75 ; MAIL1 = MANU 'QUA4' A1 A3 A5 A7 ; MAIL1 = MAIL1 ET (MANU 'TRI3' PCEN A1 A7 ); MAIL1 = MAIL1 ET (MANU 'TRI3' A1 A2 A3 ); MAIL1 = MAIL1 ET (MANU 'TRI3' A3 A4 A5 ); MAIL1 = MAIL1 ET (MANU 'TRI3' A5 A6 A7 ); MAIL2 = MAIL1; REPE BL1 3 ; MAIL2 = MAIL2 TOUR 90 PCEN ; MAIL1 = MAIL1 ET MAIL2 ; FIN BL1 ; ELIM MAIL1 0.001 ; VECX = 1.0 0.0 ; VECY = 0.0 1.0 ; NELEM = 1 ; NELEMTOT = (2 * NELEM) + 1 ; MAILTOT = MAIL1 'COUL' 'BLEU' ; IX = (-1 * NELEM) - 1 ; 'REPE' BL1 NELEMTOT ; IX = IX + 1; IY = (-1 * NELEM) - 1 ; 'REPE' BL2 NELEMTOT ; IY = IY + 1; DVECX = IX '*' VECX ; DVECY = IY '*' VECY ; DVEC = 'PLUS' DVECX DVECY ; SI ((IX NEG 0) OU (IY NEG 0)) ; MAIL2 = MAIL1 PLUS DVEC ; MAILTOT = MAILTOT 'ET' MAIL2 ; FINSI ; 'FIN' BL2 ; 'FIN' BL1 ; 'ELIM' MAILTOT 0.01 ; 'ELIM' 0.01 MAIL1 MAILTOT ; MDNS = 'EULER' ; $DOMTOT = 'MODE' MAILTOT MDNS ; $DOM1 = 'MODE' MAIL1 MDNS ; TAB1 = 'DOMA' $DOMTOT 'VF' ; TAB2 = 'DOMA' $DOM1 'VF' ; MDOMTOT = TAB1 . 'QUAF' ; MDOM1 = TAB2 . 'QUAF' ; 'ELIM' (MDOMTOT 'ET' MDOM1) 1.E-5 ; 'SI' GRAPH ; 'TRACER' ('DOMA' $DOMTOT 'MAILLAGE') 'TITRE' ('CHAINE' 'Maillage '); 'FINSI' ; *********************** **** LA TABLE PGAZ **** *********************** PGAZ = 'TABLE' ; * **** Ordre des polynoms * PGAZ . 'NORD' = 4 ; * **** Especes qui sont dans les equations d'Euler * PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ; * **** Espece qui n'y est pas * PGAZ . 'ESPNEULE' = 'N2 '; * PGAZ . 'H2 ' = 'TABLE' ; PGAZ . 'H2O ' = 'TABLE' ; PGAZ . 'N2 ' = 'TABLE' ; PGAZ . 'O2 ' = 'TABLE' ; * **** R (J/Kg/K) * PGAZ . 'H2 ' . 'R' = 4130.0 ; PGAZ . 'H2O ' . 'R' = 461.4 ; PGAZ . 'N2 ' . 'R' = 296.8 ; PGAZ . 'O2 ' . 'R' = 259.8 ; * **** Regressions polynomials * PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836 -2.37281455E-07 1.84701105E-11 ; PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05 -1.82753232E-08 2.44485692E-12 ; PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05 8.78233606E-09 -3.05514485E-13 ; PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865 2.33636971E-08 -1.53304905E-12; * **** "Enthalpies" (ou energies) de formations a OK (J/Kg) * PGAZ . 'H2 ' . 'H0K' = -4.195D6 ; PGAZ . 'H2O ' . 'H0K' = -1.395D7 ; PGAZ . 'N2 ' . 'H0K' = -2.953D5 ; PGAZ . 'O2 ' . 'H0K' = -2.634D5 ; * Scalaires passives PGAZ . 'SCALPASS' = 'MOTS' 'TAU ' 'EPS ' ; * *** Fin PGAZ * * **** ETAT CENTRE * rog = 1.0 ; ung = 0.0 ; utg = 0.0 ; retg = .4291145555695540D+04 ; yh2g = 0.01 ; yo2g = 0.2 ; yh2og = 0.3 ; taug = 1.0 ; epsg = 2.0 ; rouxg = ung '*' rog ; rouyg = utg '*' rog ; ROC1 = KCHT $DOM1 'SCAL' 'CENTRE' rog ; ROVC1 = KCHT $DOM1 'VECT' 'CENTRE' (rouxg rouyg ); RETC1 = KCHT $DOM1 'SCAL' 'CENTRE' retg ; RYH2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2 ' (rog*yh2g) ; RYO2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'O2 ' (rog*yo2g) ; RYH2OC1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rog*yh2og) ; RSCAC1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' ('MOTS' 'TAU' 'EPS') ((rog*taug) (rog*epsg)) ; * *** Le reste * rod = 1.250D-1 ; und = 0.0D0 ; utd = 0.0D0 ; retd = .3598345082089522D+01 ; yh2d = 0.2 ; yo2d = 0.4 ; yh2od = 0.1 ; taud = 11.0 ; epsd = 21.01 ; rouxd = und '*' rod ; rouyd = utd '*' rod ; ROC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' rod ; ROVC2 = KCHT $DOMTOT 'VECT' 'CENTRE' (rouxd rouyd ); RETC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' retd ; RYH2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2 ' (rod * yh2d) ; RYO2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2 ' (rod * yo2d) ; RYH2OC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rod * yh2od) ; RSCAC2 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' ('MOTS' 'TAU' 'EPS') ((rod*taud) (rod*epsd)) ; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ROC2 ROC1 ; GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' ROVC2 ROVC1 ; RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' RETC2 RETC1 ; RSCAN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' ('MOTS' 'TAU' 'EPS') RSCAC2 RSCAC1 ; RYH2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2 ' RYH2C2 RYH2C1 ; RYO2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2 ' RYO2C2 RYO2C1 ; RYH2O = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' RYH2OC2 RYH2OC1 ; RYN = RYH2 'ET' RYO2 'ET' RYH2O ; ******************************************************** *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM *** ******************************************************** MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ; * *** GRAPHIQUE DES C.I. * 'SI' GRAPH ; * 'SI' FAUX ; * *** CREATION DE CHAMELEM * CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ; CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ; CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ; CHM_RYN = 'KCHA' $DOMTOT 'CHAM' RYN ; CHM_RSN = 'KCHA' $DOMTOT 'CHAM' RSCAN ; TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RN at t=' 0.0); TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'GN at t=' 0.0); TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' 0.0); TRAC CHM_RYN MOD1 'TITR' ('CHAINE' 'RYN at t=' 0.0); TRAC CHM_RSN MOD1 'TITR' ('CHAINE' 'RSCAN at t=' 0.0); 'FINSI' ; RN0 = 'COPIER' RN ; GN0 = 'COPIER' GN ; RETN0 = 'COPIER' RETN ; RYN0 = 'COPIER' RYN ; RSCAN0 = 'COPIER' RSCAN ; * **** Parameter for the time loop * * Names of the residuum components LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' 'TAU' 'EPS' ; NOMRN = 'MOTS' 'SCAL' ; NOMVN = 'MOTS' 'UX' 'UY' ; NOMY = 'MOTS' 'H2' 'O2' 'H2O' ; NOMT = 'MOTS' 'TAU' 'EPS' ; * **** Geometrical coefficient to compute gradients * GRADR CACCA COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE' NOMRN RN ; GRADV CACCA COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE' NOMVN GN ; TPS = 0.0 ; * **** Temporal loop * 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'Methode = ' METO) ; 'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ; 'MESSAGE' ; VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSCAN ; TNM1 = 'COPIER' TN ; 'TEMPS' 'ZERO' ; 'REPETER' BL1 NITER ; * **** Primitive variables * VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSCAN TNM1 ; TNM1 = 'COPIER' TN ; 'SI' LOGSO ; GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' NOMRN RN 'GRADGEO' COEFSCAL ; GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' NOMRN PN 'GRADGEO' COEFSCAL ; GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR' NOMVN VN 'GRADGEO' COEFVECT ; GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' NOMY YN 'GRADGEO' COEFSCAL ; GRADS ALS = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' NOMT SN 'GRADGEO' COEFSCAL ; ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 2 1 $DOMTOT PGAZ RN GRADR ALR VN GRADV ALV PN GRADP ALP YN GRADY ALY SN GRADS ALS ; 'SINON' ; ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 1 1 $DOMTOT PGAZ RN VN PN YN SN ; 'FINSI' ; RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO $DOMTOT PGAZ LISTINCO ROF VITF PF YF SF ; DT_CON = SAFFAC '*' DELTAT ; * **** The time step linked to tps * DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ; * **** Total time step * DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ; * **** Increment of the variables (convection) * RESIDU = DTMIN '*' RESIDU ; DRN = 'EXCO' 'RN' RESIDU 'SCAL' ; DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ; DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ; DRYN = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU ('MOTS' 'H2' 'O2' 'H2O') ; DRSN = 'EXCO' ('MOTS' 'TAU' 'EPS') RESIDU ('MOTS' 'TAU' 'EPS') ; TPS = TPS '+' DTMIN ; RN = RN '+' DRN ; GN = GN '+' DGN ; RETN = RETN '+' DRETN ; RYN = RYN '+' DRYN ; RSCAN = RSCAN '+' DRSN ; 'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ; 'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ; 'FINSI' ; 'SI' (TPS > TFINAL) ; 'QUITTER' BL1 ; 'FINSI' ; 'FIN' BL1 ; 'TEMPS' ; * **** Les variables primitives * VN PN TN YN SN GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSCAN ; CN = (GAMN * (PN / RN)) '**' 0.5 ; * *** GRAPHIQUE DES SOLUTIONS * 'SI' GRAPH ; * *** CREATION DE CHAMELEM * CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ; CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ; CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ; CHM_RYN = 'KCHA' $DOMTOT 'CHAM' RYN ; TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' TFIN); TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'GN at t=' TFIN); TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' TFIN); TRAC CHM_RYN MOD1 'TITR' ('CHAINE' 'RYN at t=' TFIN); 'FINSI' ; * Tests: invariance par rotation RN1 = COPIER RN ; GN1 = COPI GN ; RETN1 = COPIER RETN ; RYN1 = COPIER RYN ; RSCAN1 = 'COPIER' RSCAN ; MAIL0 = EXTR RN 'MAILLAGE' ; REPETER BL1 3 ; RN1 = RN1 TOUR 90 PCEN ; GN1 = GN1 TOUR 90 PCEN ; RETN1 = RETN1 TOUR 90 PCEN ; RYN1 = RYN1 TOUR 90 PCEN ; RSCAN1 = RSCAN1 TOUR 90 PCEN ; MAIL1 = EXTR RN1 'MAILLAGE' ; ELIM 0.001 MAIL0 MAIL1 ; ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ; SI (ERR1 > 1.0D-10); ERRE 5; FINSI ; MAIL1 = EXTR GN1 'MAILLAGE' ; ELIM 0.001 MAIL0 MAIL1 ; ERR1 = MAXI ((GN1 - GN) / (RN * CN)) 'ABS' ; SI (ERR1 > 1.0D-10); ERRE 5; FINSI ; MAIL1 = EXTR RETN1 'MAILLAGE' ; ELIM 0.001 MAIL0 MAIL1 ; ERR1 = MAXI ((RETN1 - RETN) / RETN) 'ABS' ; SI (ERR1 > 1.0D-10); ERRE 5; FINSI ; MAIL1 = EXTR RYN1 'MAILLAGE' ; ELIM 0.001 MAIL0 MAIL1 ; ERR1 = MAXI ((RYN1 - RYN) / RN) 'ABS' ; SI (ERR1 > 1.0D-10); ERRE 5; FINSI ; MAIL1 = EXTR RSCAN1 'MAILLAGE' ; ELIM 0.001 MAIL0 MAIL1 ; ERR1 = MAXI ((RSCAN1 - RSCAN) / RN) 'ABS' ; SI (ERR1 > 1.0D-10); ERRE 5; FINSI ; FIN BL1 ; * Symetrie MAIL1 RN1 = MAIL0 RN SYME 'DROI' (0 0.5) (1.0 0.5) ; ELIM 0.001 MAIL0 MAIL1 ; ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ; SI (ERR1 > 1.0D-10); ERRE 5; FINSI ; FIN ;