* fichier : konv_ther_cons2.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ *********************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la **** **** solution des **** **** Equations d'Euler pour un gaz parfait **** **** OPERATEURS PRET, KONV **** **** **** **** Cas gaz multi-especes "thermally perfect" **** **** Consistence **** **** Splitting des scalaires passifs **** **** **** **** Methodes: VLH, Colella **** **** **** **** A. BECCANTINI DRN/DMT/SEMT/LTMF FEVRIER 2000 **** *********************************************************** 'OPTION' 'DIME' 2 ; 'OPTION' 'ELEM' QUA4 ; 'OPTION' 'ECHO' 0 ; 'OPTION' 'TRAC' 'X' ; * *** GRAPH * GRAPH = FAUX ; * GRAPH = VRAI ; *********************** **** LA TABLE PGAZ **** *********************** PGAZ = 'TABLE' ; PGAZ2 = 'TABLE' ; * **** Ordre des polynoms cv_i * PGAZ . 'NORD' = 1 ; PGAZ2 . 'NORD' = 1 ; * **** Especes qui sont dans les equations d'Euler * PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ; PGAZ2 . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ; * **** Espece qui n'y est pas * PGAZ . 'ESPNEULE' = 'N2 '; PGAZ2 . 'ESPNEULE' = 'N2 '; * PGAZ . 'H2 ' = 'TABLE' ; PGAZ . 'H2O ' = 'TABLE' ; PGAZ . 'N2 ' = 'TABLE' ; PGAZ . 'O2 ' = 'TABLE' ; PGAZ2 . 'H2 ' = 'TABLE' ; PGAZ2 . 'H2O ' = 'TABLE' ; PGAZ2 . 'N2 ' = 'TABLE' ; PGAZ2 . '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 ; PGAZ2 . 'H2 ' . 'R' = 4130.0 ; PGAZ2 . 'H2O ' . 'R' = 461.4 ; PGAZ2 . 'N2 ' . 'R' = 296.8 ; PGAZ2 . 'O2 ' . 'R' = 259.8 ; * **** Regressions polynomials * PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 ; PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 ; PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 ; PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 ; PGAZ2 . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 ; PGAZ2 . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 ; PGAZ2 . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 ; PGAZ2 . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 ; * **** "Enthalpies" (ou energies) de formations a OK (J/Kg) * Note: ce sont des entites numeriques PGAZ . 'H2 ' . 'H0K' = -4.195D6 ; PGAZ . 'H2O ' . 'H0K' = -1.395D7 ; PGAZ . 'N2 ' . 'H0K' = -2.953D5 ; PGAZ . 'O2 ' . 'H0K' = -2.634D5 ; PGAZ2 . 'H2 ' . 'H0K' = -4.195D6 ; PGAZ2 . 'H2O ' . 'H0K' = -1.395D7 ; PGAZ2 . 'N2 ' . 'H0K' = -2.953D5 ; PGAZ2 . 'O2 ' . 'H0K' = -2.634D5 ; * PGAZ2: transport de scalaires passifs PGAZ2 . 'SCALPASS' = 'MOTS' 'TUTU' 'TATA' ; *************************** ***** DOMAINE SPATIAL **** *************************** A1 = 0.0D0 0.0D0; A2 = 1.0D0 0.0D0; A3 = 2.0D0 0.0D0; A4 = 2.0D0 1.0D0; A5 = 1.0D0 1.0D0; A6 = 0.0D0 1.0D0; L12 = A1 'DROIT' 1 A2; L23 = A2 'DROIT' 1 A3; L34 = A3 'DROIT' 1 A4; L45 = A4 'DROIT' 1 A5; L56 = A5 'DROIT' 1 A6; L61 = A6 'DROIT' 1 A1; L25 = A2 'DROIT' 1 A5; DOM10 = 'DALL' L12 L25 L56 L61 'PLANE'; DOM20 = 'DALL' L23 L34 L45 ('INVERSE' L25) 'PLANE'; * *** Point ou on controlle la consistence * P10 = 1.0 0.5; * *** Etats gauche et droite * rogd = 1.412 ; pgd = 1.17D6 ; ungd = 8.34 ; utgd = 18.2 ; yh2gd = 0.1 ; yh2Ogd = 0.2 ; yN2gd = 0.15 ; yO2gd = 1 '-' (yh2gd '+' yh2Ogd '+' yn2gd) ; tutugd = 14.9 ; tatagd = 23.8 ; * **** tg * tgd = pgd '/' (rogd * ((yh2gd * (PGAZ . 'H2 ' . 'R')) '+' (yh2Ogd * (PGAZ . 'H2O ' . 'R')) '+' (yn2gd * (PGAZ . 'N2 ' . 'R')) '+' (yo2gd * (PGAZ . 'O2 ' . 'R')))) ; * **** gamg, gamd , htg , htd * * A0H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 1 ; A1H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 2 ; A0H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 1 ; A1H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 2 ; A0N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 1 ; A1N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 2 ; A0O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 1 ; A1O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 2 ; T = tgd ; T2 = T * T ; ETGD = (yh2gd * ((A0H2 * T) '+' (A1H2 * T2 '/' 2))); ETGD = ETGD '+' (yH2Ogd * ((A0H2O * T) '+' (A1H2O * T2 '/' 2))); ETGD = ETGD '+' (yN2gd * ((A0N2 * T) '+' (A1N2 * T2 '/' 2))); ETGD = ETGD '+' (yO2gd * ((A0O2 * T) '+' (A1O2 * T2 '/' 2))) ; ecingd = 0.5D0 '*' ((ungd '*' ungd) '+' (utgd '*' utgd)); retgd = rogd '*' (ETGd '+' ecingd) ; * *** flux en (n,t) * f1gd = ungd '*' rogd ; f2gd = (rogd '*' (ungd '*' ungd)) '+' pgd ; f3gd = rogd '*' ungd '*' utgd ; f4gd = ungd '*' (retgd '+' pgd); f5gd = ungd '*' rogd '*' yh2gd ; f6gd = ungd '*' rogd '*' yo2gd ; f7gd = ungd '*' rogd '*' yh2ogd ; ftutugd = ungd '*' rogd '*' tutugd ; ftatagd = ungd '*' rogd '*' tatagd ; **************************************************** **************************************************** ******** Boucle sur les angles ********* **************************************************** **************************************************** DANGLE = 360 '/' 7.15; ANGLE = 11.3 ; 'REPETER' BLMET 2 ; SI (&BLMET 'EGA' 1) ; METO = 'VLH' ; 'FINSI' ; SI (&BLMET 'EGA' 2) ; METO = 'SS' ; 'FINSI' ; 'REPETER' BLOC 8; * *** Rotation * ANGLE = ANGLE '+' DANGLE; ORIG = 0.0D0 0.0D0; MESSAGE; MESSAGE (CHAIN 'Angle de rotation= ' ANGLE); MESSAGE; DOM1 = DOM10 'TOURNER' ANGLE ORIG; DOM2 = DOM20 'TOURNER' ANGLE ORIG; P1 = P10 'TOURNER' ANGLE ORIG; 'ELIMINATION' (DOM1 ET DOM2) 1D-6; DOMTOT = DOM1 ET DOM2; $DOMTOT = 'MODELISER' DOMTOT 'EULER'; $DOM1 = 'MODELISER' DOM1 'EULER'; $DOM2 = 'MODELISER' DOM2 'EULER'; TDOMTOT = 'DOMA' $DOMTOT 'VF'; TDOM1 = 'DOMA' $DOM1 'VF'; TDOM2 = 'DOMA' $DOM2 'VF'; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMTOT = TDOMTOT . 'QUAF' ; 'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) 1.E-6 ; 'SI' GRAPH; 'TRACER' (('DOMA' $DOMTOT 'MAILLAGE' ) 'ET' ('DOMA' $DOMTOT 'FACEL') 'ET' P1) 'TITRE' 'Domaine et FACEL'; 'FINSI' ; * **** Redefinition de P1 dans $DOMTOT 'FACE' * P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1; *********************** **** Les CHPOINTs **** *********************** uxgd = (ungd '*' ('COS' ANGLE)) '-' (utgd '*' ('SIN' ANGLE)); uygd = (ungd '*' ('SIN' ANGLE)) '+' (utgd '*' ('COS' ANGLE)); RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rogd; RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rogd;; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2); VN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (uxgd uygd); VN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (uxgd uygd); VN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VN1 'ET' VN2); PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pgd; PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pgd; PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2); * attention: PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ; YH2 = ('KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2 ' yh2gd) 'ET' ('KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'H2 ' yh2gd) ; YO2 = ('KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'O2 ' yo2gd) 'ET' ('KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'O2 ' yo2gd) ; YH2O = ('KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2O ' yh2ogd) 'ET' ('KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'H2O ' yh2ogd) ; YN = YH2 'ET' YO2 'ET' YH2O ; * SCAN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ; SCAN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' (PGAZ2 . 'SCALPASS') (tutugd tatagd) ; SCAN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' (PGAZ2 . 'SCALPASS') (SCAN1 'ET' SCAN2); *************************** **** L'operateur PRET**** *************************** ORDESP = 1; ORDTEM = 1; ROF VITF PF YF = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN VN PN YN ; ROF2 VITF2 PF2 YF2 SCAF = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN VN PN YN SCAN ; GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1; GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1; REFGEOP1 = 'REDU' VITF GEOP2; * **** Orientation de la normal n de castem par raport a la * notre; t est par consequence * NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1; NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1; 'SI' (('ABS' NCOS) > ('ABS' NSIN)); ORIENT = ('COS' ANGLE) '/' NCOS; 'SINON'; ORIENT = ('SIN' ANGLE) '/' NSIN; 'FINSI' ; ORIENT = 'SIGN' ORIENT; ************************** ************************** **** 'KONV', TEST VLH **** ************************** ************************** LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ; TFLUX DELTAT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' METO $DOMTOT PGAZ LISTINCO ROF VITF PF YF ; * **** Les flux aux interfaces sont dans de * CHPOINT FACE * FLUX1 = 'EXCO' 'RN' TFLUX ; FLUX2X = 'EXCO' 'RUX' TFLUX ; FLUX2Y = 'EXCO' 'RUY' TFLUX ; FLUX3 = 'EXCO' 'RETN' TFLUX ; FLUX4 = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') TFLUX ('MOTS' 'H2' 'O2' 'H2O') ; FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE)); FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE)); * **** Test 'KONV' * f1 = 'EXTRAIRE' FLUX1 'SCAL' P1; ERRO = 'ABS' (1D-12 '*' f1gd); LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO; f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1; ERRO = 'ABS' (1D-12 '*' f2gd); LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO; LOGI1 = LOGI1 'ET' LOGI2; f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1; ERRO = 'ABS' (1D-11 '*' f3gd) ; LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO; LOGI1 = LOGI1 'ET' LOGI2; f4 = 'EXTRAIRE' FLUX3 'SCAL' P1; ERRO = 'ABS' (1D-12 '*' f4gd) ; LOGI2 = ('ABS' ((f4 '*' ORIENT) '-' f4gd)) < ERRO; LOGI1 = LOGI1 'ET' LOGI2; f5 = 'EXTRAIRE' FLUX4 'H2 ' P1; ERRO = 'ABS' (1D-12 '*' f5gd) ; LOGI2 = ('ABS' ((f5 '*' ORIENT) '-' f5gd)) < ERRO; LOGI1 = LOGI1 'ET' LOGI2; f6 = 'EXTRAIRE' FLUX4 'O2 ' P1; ERRO = 'ABS' (1D-12 '*' f6gd) ; LOGI2 = ('ABS' ((f6 '*' ORIENT) '-' f6gd)) < ERRO; LOGI1 = LOGI1 'ET' LOGI2; f7 = 'EXTRAIRE' FLUX4 'H2O ' P1; ERRO = 'ABS' (1D-12 '*' f7gd) ; LOGI2 = ('ABS' ((f7 '*' ORIENT) '-' f7gd)) < ERRO; LOGI1 = LOGI1 'ET' LOGI2; 'SI' ('NON' LOGI1); 'MESSAGE' ; 'MESSAGE' 'OPERATEUR KONV'; 'MESSAGE' ('CHAINE' METO); 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'df1 ' ('ABS' (f1gd '-' (ORIENT '*'f1))) ' f1 ' f1); 'MESSAGE' ('CHAINE' 'df2 ' ('ABS' (f2gd '-' (ORIENT '*'f2))) ' f2 ' f2); 'MESSAGE' ('CHAINE' 'df3 ' ('ABS' (f3gd '-' (ORIENT '*'f3))) ' f3 ' f3); 'MESSAGE' ('CHAINE' 'df4 ' ('ABS' (f4gd '-' (ORIENT '*'f4))) ' f4 ' f4); 'MESSAGE' ('CHAINE' 'df5 ' ('ABS' (f5gd '-' (ORIENT '*'f5))) ' f5 ' f5); 'MESSAGE' ('CHAINE' 'df6 ' ('ABS' (f6gd '-' (ORIENT '*'f6))) ' f6 ' f6); 'MESSAGE' ('CHAINE' 'df7 ' ('ABS' (f7gd '-' (ORIENT '*'f7))) ' f7 ' f7); 'ERREUR' 5; 'FINSI' ; ************************************************ ***** Splitting des scalaires passifs ********** ************************************************ * LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' 'TUTU' 'TATA' ; TFLUX2 DELTAT = 'KONV' 'VF' 'PERFTEMP' 'FLUX' METO $DOMTOT PGAZ2 LISTINCO ROF2 VITF2 PF2 YF2 SCAF ; TFLUX3 = 'EXCO' ( 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ) TFLUX2 ( 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ) ; ERRO = 'MAXIMUM' (TFLUX3 '-' TFLUX) 'ABS' ; 'SI' (ERRO > 1.0D-15); 'MESSAGE' 'Probleme dans les transport de scalaires passifs' ; 'ERREUR' 5; 'FINSI' ; * **** Test 'KONV' * fceltutu = 'EXTRAIRE' TFLUX2 'TUTU' P1; LOGI1 = ('ABS' ((fceltutu '*' ORIENT) '-' ftutugd)) < ('ABS' (1.0D-12 * fceltutu)) ; fceltata = 'EXTRAIRE' (TFLUX2) 'TATA' P1; LOGI2 = ('ABS' ((fceltata '*' ORIENT) '-' ftatagd)) < ('ABS' (1.0D-12 * fceltata)) ; LOGI1 = LOGI1 'ET' LOGI2; 'SI' ('NON' LOGI1); 'MESSAGE' ; 'MESSAGE' 'OPERATEUR KONV'; 'MESSAGE' 'Transport de scalaires passifs' ; 'ERREUR' 5; 'FINSI' ; **************************************************** **************************************************** ******** Fin boucle sur les angles ********* **************************************************** **************************************************** 'FIN' BLOC; 'FIN' ;