* fichier : rut_tg_2.dgibi ************************************************************************ ************************************************************************ ************************************************************************ **************** CAS-TEST RUT ************************************* ************************************************************************ * Date de creation : 16 decembre 2002 * Auteur : Alexandre Bleyer * Commentaire : ce cas-test est de coupe en plusieurs * parties : procedures, definition du maillage, de la table PGAZ * et definition des conditons initiales ************************************************************************ * Liste des procedures : CONS, DETO, EXEX, PERSO, PREPADET, VALP ************************************************************************ *$$$$ CONS *************************************************************** ***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES ***** *************************************************************** 'DEBPROC' CONS RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' YN*'CHPOINT' $DM*'MMODEL' PGAZ*'TABLE' ; * NOMCOM = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE'); YTOT = YTOT '+' YCEL ; RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') * YCEL) ; A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; 'SI' ('EGA' &BL1 1) ; 'SINON' ; RYN = RYN 'ET' 'FINSI' ; 'FIN' BL1 ; NOMCOM = PGAZ . 'ESPNEULE' ; A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL) ; RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') '*' YCEL); TN = PN '/' (RN '*' RTOT) ; T2 = TN '*' TN ; T3 = T2 '*' TN ; T4 = T3 '*' TN ; T5 = T4 '*' TN ; ETHER = (A0 * TN) '+' ((A1 '/' 2.0) * T2) '+' ((A2 '/' 3.0) * T3) '+' ((A3 '/' 4.0) * T4) '+' ((A4 '/' 5.0) * T5) ; REN = (RN '*' ETHER) '+' ECIN ; 'FINPROC' GN REN RYN ; *$$$$ DETO 'ARGUMENT' PGAZ*'TABLE' TN*'CHPOINT' RYN*'CHPOINT' DTI*'FLOTTANT' ; ************************ **** PROCEDURE DETO **** ************************ * * dRYO2/dT = MO2 * A * (T ** -B) EXP( - TA / T) * (RYH2 ** C) * (RYO2) * **** Les parametres physiques * * Ts = temperature de seuil * TA = 8310.; A = 1.1725D14; B = 0.91; CON = 2; H0H2 = PGAZ . 'H2 ' . 'H0K' ; H0O2 = PGAZ . 'O2 ' . 'H0K' ; H0H2O = PGAZ . 'H2O ' . 'H0K' ; RYH2 RYO2 RYH2O DELTAE = 'FLAM' 'ARRHENIU' TS A B CON TA H0H2 H0O2 H0H2O DTI RYH2 RYO2 RYH2O TN ; RY = RYH2 'ET' RYO2 'ET' RYH2O ; 'FINPROC' RY DELTAE ; **************************** **** FIN PROCEDURE DETO **** **************************** *$$$$ EXEX 'DEBPROC' EXEX RV*'TABLE' ; **************************************************************** **************************************************************** ***** Parameters for the computations ****** **************************************************************** **************************************************************** * Upwind scheme METO = RV.'CI'.'METO' ; * Iterations * Final time * Safety Factor for the time step * Second order reconstruction? TPS = RV . 'PASDETPS' . 'TINI' ; 'SI' (NTPS '>' 1) ; 'FINSI'; NITER = RV . 'PASDETPS' . 'NITMA' ; TFINAL = RV . 'PASDETPS' . 'TFINA' ; SAFFAC = RV . 'PASDETPS' . 'SAFFAC' ; LOGSO = 'EGA' (RV.'ORDREESP') 2 ; GEO = RV . 'GEO' ; INCO = RV . 'INCO' ; PGAZ = RV . 'PGAZ' ; $DOMTOT= GEO . 'MODTOT' ; * **** Chemical time step * * We compute d(RHY2)/dT at the von Neumann state * * Molar reaction rate * TVN = 2000 ; TA = 8310.; A = 1.1725D14; B = 0.91; CON = 2; H0H2 = PGAZ . 'H2 ' . 'H0K' ; H0O2 = PGAZ . 'O2 ' . 'H0K' ; H0H2O = PGAZ . 'H2O ' . 'H0K' ; RYN = INCO . 'RYN' ; GN = INCO . 'GN' ; RN = INCO . 'RN' ; RETN = INCO . 'RETN' ; IGNI1 = 2 ; MOD_I = 'CHAIN' 'MOD' IGNI1; $DOM2 = GEO . MOD_I; OMEGA = A '*' (RYH2 '**' CON) '*' RYO2 '*' (TVN '**' (-1.0 '*' B)) '*' ('EXP' (-1.0 '*' TA '*' (TVN '**' -1))) ; * OMEGA in (mol/m3/s) * OMEGAH2 in (kg/m3/s) OMEGAH2 = OMEGA '*' 2 '*' 2E-3 ; DELTATC = RYH2 '/' OMEGAH2 ; * NB DELTATC << DELTAT_CON * **** Zone d activation = etat AICC * IGNI1 = RV . 'MODETO' .'IGNIREG' ; MOD_I = 'CHAIN' 'MOD' IGNI1; $DOM1 = GEO . MOD_I; ; RYH2 RYO2 RYH2O DELTARE = 'FLAM' 'ARRHENIU' TS A B CON TA H0H2 H0O2 H0H2O (1.0D5 '*' DELTATC) RYH2 RYO2 RYH2O TN ; RYN = RYH2 'ET' RYO2 'ET' RYH2O ; RETN = RETN '+' DELTARE ; RN0 = 'COPIER' RN ; GN0 = 'COPIER' GN ; RETN0 = 'COPIER' RETN ; RYN0 = 'COPIER' RYN ; * **** Parameter for the time loop * * Names of the residuum components * **** Geometrical coefficient to compute gradients * GRADR CACCA COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE' GRADV CACCA COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE' * **** The chemistry time step is too little * We take DT_CHEM > DELTATC * * DT_CHEM = 50 '*' DELTATC ; * **** Temporal loop * 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'Methode = ' METO) ; 'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ; 'MESSAGE' ; TNM1 = 'COPIER' TN ; 'TEMPS' 'ZERO' ; 'REPETER' BL1 NITER ; * **** Primitive variables * VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN TNM1 ; TNM1 = 'COPIER' TN ; 'SI' LOGSO ; GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR' GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' $DOMTOT PGAZ RN GRADR ALR VN GRADV ALV PN GRADP ALP YN GRADY ALY ; 'SINON' ; $DOMTOT PGAZ RN VN PN YN ; 'FINSI' ; $DOMTOT PGAZ LISTINCO ROF VITF PF YF ; DT_CON = SAFFAC '*' DELTAT ; * **** The time step linked to tps * * **** Total time step * DTIMP = RV.'PASDETPS'.'DTIMP'; * **** Increment of the variables (convection) * RESIDU = DTMIN '*' RESIDU ; TPS = TPS '+' DTMIN ; RV.'PASDETPS'.'LTPS' = RV.'PASDETPS'.'LTPS' 'ET' RV.'PASDETPS'.'TPS' = TPS ; RN = RN '+' DRN ; GN = GN '+' DGN ; RETN = RETN '+' DRETN ; RYN = RYN '+' DRYN ; * **** Increment of the variables (chemical source) * VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN TNM1 ; RETN = RETN '+' DRETN ; RV.'INCO'.'PN' = PN ; RV.'INCO'.'TN' = TN ; RV.'INCO'.'RN' = RN ; RV.'INCO'.'VN' = VN ; RV.'INCO'.'YN' = YN ; RV.'INCO'.'GAMMAN' = GAMMAN ; RV.'INCO'.'GN' = GN ; RV.'INCO'.'RETN' = RETN ; RV.'INCO'.'RYN' = RYN ; 'SI' ('EXISTE' RV 'TABC' ); RV . 'TABC' = TABC ; 'FINSI' ; 'SI' ('EXISTE' RV 'PRCPERSO' ); FREQ0 = RV.'PRCPERSO'.'ARG3' ; 'SI' (((&BL1 '/' FREQ0) '*' FREQ0) 'EGA' &BL1) ; (RV.'PRCPERSO'.'ARG2') ; 'FINSI' ; 'FINSI' ; 'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ; 'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ; 'FINSI' ; 'SI' (TPS '>EG' TFINAL) ; 'QUITTER' BL1 ; 'FINSI' ; 'FIN' BL1 ; 'TEMPS' ; RV.'PASDETPS'.'TINI' = TPS ; 'FINPROC' RV; *$$$$ PERSO 'DEBP' PERSO RV*'TABLE' $DOMTRA*'MMODEL' TRACE*'ENTIER'; INCO = RV . 'INCO' ; GEO = RV . 'GEO' ; TPS = RV . 'PASDETPS' . 'TPS' ; * X0 Y0 = 'COORD' DOMT1 ; DTX = (Xmax - Xmin)/5. ; DTY = (Ymax - Ymin)/3. ; Pt1 = (Xmin - DTX) (Ymin - DTY) ; Pt2 = (Xmax + DTX) (Ymin - DTY) ; Pt3 = (Xmax + DTX) (Ymax + DTY) ; Pt4 = (Xmin - DTX) (Ymax + DTY) ; TRACHA = FAUX ; TRACAP = FAUX ; 'SI' (TRACE 'EGA' 1); TRACHA = VRAI ; 'FINSI'; 'SI' (TRACE 'EGA' 2); TRACAP = VRAI ; 'FINSI'; 'SI' (TRACE 'EGA' 3); TRACHA = VRAI ; TRACAP = VRAI ; 'FINSI'; 'SI' TRACHA ; TIT1 = 'CHAIN' 'PRESSION' TPS ; * TIT1 = 'CHAIN' 'TEMPERATURE' TPS ; * *AB = KCHA $DOMT1 INCO.'RN' 'CHAM' ; *TIT1 = 'CHAIN' 'DENSITE' TPS ; *TRAC AB $DOMT1 cont1 'TITR' TIT1 ; * *YN = INCO.'YN' ; *AA = 'EXCO' 'H2 ' YN ; *AB = KCHA $DOMT1 AA 'CHAM' ; *TIT1 = 'CHAIN' 'YH2' TPS ; *TRAC AB $DOMT1 cont1 'TITR' TIT1 ; * *AA = 'EXCO' 'O2 ' YN ; *AB = KCHA $DOMT1 AA 'CHAM' ; *TIT1 = 'CHAIN' 'YO2' TPS ; *TRAC AB $DOMT1 cont1 'TITR' TIT1 ; * *AA = 'EXCO' 'H2O ' YN ; *AB = KCHA $DOMT1 AA 'CHAM' ; *TIT1 = 'CHAIN' 'YH2O' TPS ; *TRAC AB $DOMT1 cont1 'TITR' TIT1 ; * TIT1 = 'CHAIN' 'VITESSE' TPS ; * 'FINSI'; 'SI' TRACAP ; TABC = RV . 'TABC' ; * *- Initialisations * LEGEND = 'TABLE' ; TEVOL = 'TABLE' ; 'REPETER' BLO10 NBCOMP ; IPOS2 = 1 ; IPOS0 = 0 ; I100 = 0 ; 'REPETER' BLO100 NBEVOL ; I100 = I100 + 1 ; 'SINON' ; SUPP1 = 'CHAINE' 'CAPT' I100; 'FINSI'; IPOS1 = I100 - (7 * (I100 / 7)) ; 'SI' (IPOS1 EGA 0) ; IPOS0 = 1 ; IPOS1 = 7 ; 'FINSI' ; IPOS3 = (IPOS1 + 1) - (2 * ((IPOS1 + 1) / 2)) ; 'SI' (IPOS3 EGA 0) ; IPOS2 = IPOS2 + 1 ; 'SI' (IPOS2 EGA 5) ; IPOS2 = 1 ; IPOS3 = IPOS1 ; 'FINSI' ; 'FINSI' ; LEGEN1 = 'CHAIN' 'MARQ ' IMARQ1 ; LEGEN2 = 'CHAIN' SUPP1 ; * 'SI' (IPOS1 EGA 1) ; TEVOL . LCOMP1 = EVOL1 ; LEGEND . LCOMP1 = TABLE ; LEGEND . LCOMP1 . 'TITRE' = TABLE ; 'SINON' ; TEVOL . LCOMP1 = (TEVOL . LCOMP1) ET EVOL1 ; 'FINSI' ; 'SI' ((IPOS0 EGA 1) 'OU' (I100 'EGA' NBEVOL)) ; IPOS0 = 0 ; 'FINSI' ; * 'FIN' BLO100 ; 'FIN' BLO10 ; 'FINSI' ; * FINP ; *** Debut proc perso *$$$$ PREPADET 'DEBPROC' PREPADET CONDINI*'TABLE' PGAZ*'TABLE' TABC/'TABLE'; ERROR = 0 ; *** DEFINITION DU PARAMETRE EPSI POUR ELIM DU MAILLAGE 'SINON'; 'FINSI'; *** DEFINITION DU PARAMETRE EPSI POUR ELIM DU MAILLAGE 'SI' ('NON' ('EXISTE' CONDINI 'GRAPH')); CONDINI . 'GRAPH' = FAUX ; 'FINSI'; *** INITILISATION DE TABLES TEMPORAIRES OU NON INI = 'TABLE' ; INCO = 'TABLE' 'INCO' ; GEO = 'TABLE' 'GEO'; FRACY = 'TABLE' ; *** CREATION DES OBJETS MODELES *** ils sont sauvegardes dans la table GEO de RV MAILTOT = CONDINI . 'MAILTOT' ; *** Parametres verifiant si l utilisateur a defini *** - soit des fractions massiques *** - soit des fractions volumiques *** Pour le calcul on n utilise uniquement des fractions *** massiques. *** Si l utilisateur a defini des fractions volumiques *** il faut se ramener a des fractions massiques (pre-calculs) LOGX = FAUX ; LOGY = FAUX ; 'REPETER' BCL1 CONDINI . 'NBDOM'; MODE_I = 'CHAIN' 'MOD' &BCL1; MQUA_I = 'CHAIN' 'MQUA' &BCL1; DOM_I = 'CHAIN' 'DOM' &BCL1; M_I = 'CHAIN' 'MAIL' &BCL1; MAIL_I = CONDINI.DOM_I.'MAIL' ; GEO.MODE_I = $MAIL_I ; GEO.MQUA_I = MMAIL_I ; 'SI' (&BCL1 'EGA' 1); LOGX = 'EXISTE' CONDINI.DOM_I 'XH2' ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XO2') ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XN2') ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XH2O'); LOGY = 'EXISTE' CONDINI.DOM_I 'YH2' ; LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YO2') ; LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YN2') ; LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YH2O'); 'SINON' ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XH2') ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XO2') ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XN2') ; LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XH2O'); LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YH2') ; LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YO2') ; LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YN2') ; LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YH2O'); 'FINSI'; 'FIN' BCL1; GEO.'MODTOT' = $MAILTOT ; GEO.'MQUATOT' = MMAILTOT ; *** Verification du bon remplissage de la table de donnees 'SI' (LOGX 'ET' ('NON' LOGY)); 'MESS' 'C est bien ! Vous n avez defini que des fractions volumiques'; 'SINON'; 'SI' (LOGY 'ET' ('NON' LOGX)); 'MESS' 'C est bien ! Vous n avez defini que des fractions massiques'; 'SINON'; 'MESS' 'Vous ne devez definir uniquement des fractions massiques'; 'MESS' 'En aucun cas, un melange des deux !!!'; ERROR = ERROR + 1 ; 'FINSI'; 'FINSI'; 'SI' LOGX ; FRACY.'MH2' = 2. '*' 1.00797E-3 ; FRACY.'MO2' = 2. '*' 15.9994E-3 ; FRACY.'MH2O' = (FRACY.'MH2') '+' (0.5 '*' (FRACY.'MO2')) ; FRACY.'MN2' = 2. '*' 14.0067E-3 ; 'FINSI'; 'REPETER' BCL1 CONDINI . 'NBDOM'; DOM_I = 'CHAIN' 'DOM' &BCL1; FRACY.DOM_I = 'TABLE'; 'SI' LOGX; 'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'T')); 'MESS' 'Comme vous avez defini les fractions volumiques,'; 'MESS' 'vous devez definir la temperature !'; ERROR = ERROR + 1 ; 'FINSI'; *** CALCUL DES MASSES TOTALES!!! FRACY.DOM_I.'MTOT' = ((CONDINI.DOM_I.'XH2') '*' (FRACY.'MH2')) '+' ((CONDINI.DOM_I.'XO2') '*' (FRACY.'MO2')) '+' ((CONDINI.DOM_I.'XH2O') '*' (FRACY.'MH2O')) '+' ((CONDINI.DOM_I.'XN2') '*' (FRACY.'MN2')) ; *** FIN CALCUL DES MASSES TOTALES 'FINSI'; ***************** 'SI' LOGY; 'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'RHO')); 'MESS' 'Comme vous avez defini les fractions massiques,'; 'MESS' 'vous devez definir la densite du melange !'; ERROR = ERROR + 1 ; 'FINSI'; 'FINSI'; ***************** 'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'P')); 'MESS' 'Il manque l initialisation de la pression'; ERROR = ERROR + 1 ; 'FINSI'; ***************** 'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'V')); 'MESS' 'Il manque l initialisation de la vitesse'; ERROR = ERROR + 1 ; 'FINSI'; 'FIN' BCL1; ***************** 'SI' (ERROR 'NEG' 0); 'MESS' 'Il y a des problèmes d initialisation !'; 'FINSI'; ****** TRANSFORMATION DES FRACTIONS VOLUMIQUES EN MASSIQUES 'SI' LOGX ; 'REPETER' BCL1 CONDINI . 'NBDOM'; DOM_I = 'CHAIN' 'DOM' &BCL1; FRACY.DOM_I.'YH2' = ((CONDINI.DOM_I.'XH2') '*' (FRACY.'MH2')) '/' (FRACY.DOM_I.'MTOT') ; FRACY.DOM_I.'YO2' = ((CONDINI.DOM_I.'XO2') '*' FRACY.'MO2') '/' (FRACY.DOM_I.'MTOT') ; FRACY.DOM_I.'YH2O'= ((CONDINI.DOM_I.'XH2O') '*' FRACY.'MH2O') '/' (FRACY.DOM_I.'MTOT') ; FRACY.DOM_I.'YN2' = ((CONDINI.DOM_I.'XN2') '*' FRACY.'MN2') '/' (FRACY.DOM_I.'MTOT') ; FRACY.DOM_I.'RTOT'= ((FRACY.DOM_I.'YH2') '*' (PGAZ.'H2'.'R')) '+' ((FRACY.DOM_I.'YO2') '*' (PGAZ.'O2'.'R')) '+' ((FRACY.DOM_I.'YH2O') '*' (PGAZ.'H2O'.'R')) '+' ((FRACY.DOM_I.'YN2') '*' (PGAZ.'N2'.'R')) ; FRACY.DOM_I.'RHO' = (CONDINI.DOM_I.'P') '/' ((FRACY.DOM_I.'RTOT') '*' (CONDINI.DOM_I.'T')); 'FIN' BCL1; 'FINSI'; **** INITIALISATION DES CHAMPOINTS POUR LE CALCUL 'REPETER' BCL1 CONDINI . 'NBDOM'; MODE_I = 'CHAIN' 'MOD' &BCL1; DOM_I = 'CHAIN' 'DOM' &BCL1; RHO_I = 'CHAIN' 'RHO' &BCL1; P_I = 'CHAIN' 'P' &BCL1; V_I = 'CHAIN' 'V' &BCL1; YH2_I = 'CHAIN' 'YH2' &BCL1; YO2_I = 'CHAIN' 'YO2' &BCL1; YH2O_I = 'CHAIN' 'YH2O' &BCL1; CONDINI.DOM_I.'P'; CONDINI.DOM_I.'V'; 'SI' LOGY; CONDINI.DOM_I.'RHO'; 'SINON'; 'SI' LOGX; FRACY.DOM_I.'RHO'; 'FINSI'; 'FINSI'; 'SI' (&BCL1 'EGA' 1); YH2= 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' YO2= 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' YH2O= 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' 'SINON'; 'FINSI'; 'FIN' BCL1; *** On supprime les tables qui ne nous sont plus necessaires 'OUBLIER' INI ; 'OUBLIER' FRACY ; *** YN = YH2 'ET' YO2 'ET' YH2O ; INCO . 'YN' = YN ; 'SI' CONDINI.'GRAPH'; 'SINON'; 'CACH' 'TITR' 'RO' ; 'CACH' 'TITR' 'YH2' ; 'CACH' 'TITR' 'YO2' ; 'CACH' 'TITR' 'YH2O'; 'CACH' 'TITR' 'P' ; 'FINSI'; 'FINSI'; GN REN RYN = CONS RN VN PN YN GEO.'MODTOT' PGAZ ; *************************** **** La table PASDETPS **** *************************** PASDETPS = 'TABLE' 'PASDETPS' ; 'SINON'; 'FINSI'; 'SI' ('EXISTE' CONDINI 'NITMA'); PASDETPS.'NITMA' = CONDINI . 'NITMA' ; 'SINON'; PASDETPS.'NITMA' = 10000 ; 'FINSI'; 'SI' ('EXISTE' CONDINI 'TINI'); PASDETPS.'TINI' = CONDINI . 'TINI' ; 'SINON'; PASDETPS.'TINI' = 0.0 ; 'FINSI'; PASDETPS.'TPS' = PASDETPS.'TINI' ; 'SI' ('EXISTE' CONDINI 'TFINA'); PASDETPS.'TFINA' = CONDINI . 'TFINA' ; 'SINON'; 'MESS' 'Vous devez entrer un temps final'; ERROR = ERROR + 1 ; 'FINSI'; 'SI' ('EXISTE' CONDINI 'DTIMP'); PASDETPS.'DTIMP' = CONDINI . 'DTIMP' ; 'SINON'; 'MESS' 'Vous devez entrer un pas de temps'; ERROR = ERROR + 1 ; 'FINSI'; * **** * 'SI' ('EXISTE' CONDINI 'SAFFAC'); SAFFAC = CONDINI . 'SAFFAC' ; 'SINON'; SAFFAC = 0.7 ; 'MESS' 'On prend le facteur de securite pour le pas de temps'; 'FINSI'; PASDETPS.'SAFFAC' = SAFFAC ; 'SI' ('EXISTE' CONDINI 'METO'); METO = CONDINI . 'METO' ; 'SINON'; METO = 'VLH '; 'MESS' 'Methode VLH par defaut !!!'; 'FINSI'; CONDINI.'METO' = METO ; RV = 'TABLE' ; RV . 'INCO' = INCO ; 'SI' ('EXISTE' 'TABC'); RV . 'TABC' = TABC ; 'FINSI'; RV . 'CI' = CONDINI ; RV . 'GEO' = GEO ; RV . 'PASDETPS'= PASDETPS ; RV . 'INCO' . 'RN' = 'COPIER' RN ; RV . 'INCO' . 'GN' = 'COPIER' GN ; RV . 'INCO' . 'RETN' = 'COPIER' REN ; RV . 'INCO' . 'RYN' = 'COPIER' RYN ; * *** Le gaz * RV . 'PGAZ' = PGAZ ; * **** Ordre en espace * 'SI' ('EXISTE' CONDINI 'ORDREESP'); IE = CONDINI . 'ORDREESP' ; 'SINON'; IE = 1 ; 'FINSI'; RV . 'ORDREESP' = IE ; * **** Ordre en temps * 'SI' ('EXISTE' CONDINI 'ORDRETPS'); IT = CONDINI . 'ORDRETPS' ; 'SINON'; IT = 1 ; 'FINSI'; RV . 'ORDRETPS' = IT ; * **** Definition du modele de combustion/detonation * 'SI' ('EXISTE' CONDINI 'MODETO'); MODETO = CONDINI . 'MODETO' ; 'SINON'; MODETO = 'ARRHENIU' ; 'FINSI'; RV . 'MODETO' .'NOM' = MODETO ; **** Si le modele de combustion est CREBCOM : 'SI' ('EGA' MODETO 'CREBCOM'); * numero du domaine d'ignition 'SI' ('EXISTE' CONDINI 'IGNIREG'); RV . 'MODETO' .'IGNIREG' = CONDINI . 'IGNIREG' ; 'SINON'; RV . 'MODETO' .'IGNIREG' = 1 ; 'MESS' 'On prend le 1er domaine comme domaine d ignition'; 'FINSI'; * critere Crebcom (chimie) 'SI' ('EXISTE' CONDINI 'EPS_CC'); RV . 'MODETO' .'EPS_CC' = CONDINI . 'EPS_CC' ; 'SINON'; RV . 'MODETO' .'EPS_CC' = 0.7 ; 'FINSI'; * critere Crebcom (vitesse) 'SI' ('EXISTE' CONDINI 'K0'); RV . 'MODETO' .'K0' = CONDINI . 'K0' ; 'SINON'; RV . 'MODETO' .'K0' = 500.0 ; 'FINSI'; 'FINSI'; **** Si le modele de combustion est ARRHENIUS : 'SI' ('EGA' MODETO 'ARRHENIU'); * numero du domaine d'ignition 'SI' ('EXISTE' CONDINI 'IGNIREG'); RV . 'MODETO' .'IGNIREG' = CONDINI . 'IGNIREG' ; 'SINON'; RV . 'MODETO' .'IGNIREG' = 1 ; 'MESS' 'On prend le 1er domaine comme domaine d ignition'; 'FINSI'; 'FINSI'; * **** * 'SI' (ERROR 'NEG' 0); 'MESS' 'Il y a des problèmes d initialisation !'; 'FINSI'; 'FINPROC' RV; *$$$$ VALP ***************************************************** ***************************************************** ** PROCEDURE VALP POUR CALCULER LES VALEURS DE ** ** PRESSION EN DES POINTS DONNES ET A DES TEMPS ** ** DONNES ** ***************************************************** ***************************************************** TABC = RV . 'TABC' ; MODTOT = RV . 'GEO' . 'MODTOT' ; PTS = TABC . 'LPOINTS' ; LTPS = TABC . 'LTPS' ; * On etudie uniquement la pression * peut evoluer par la suite !!! * Initialisation des tables 'SI' (NTPS EGA 0 ); 'REPETER' BOU1 M1; 'SINON' ; 'FINSI'; 'SI' (&BOU1 'EGA' 1); TABC.'RPOINTS' = PT2 ; 'SINON'; TABC.'RPOINTS' = RV.'TABC'.'RPOINTS' 'ET' PT2; 'FINSI'; 'REPETER' BOU2 N1; 'FIN' BOU2; 'FIN' BOU1; 'SI' RV.'CI'.'GRAPH' ; 'SINON'; 'FINSI'; 'FINSI'; 'FINSI'; 'REPETER' BOU1 M1; 'SINON' ; 'FINSI'; 'REPETER' BOU2 N1; *** captage de la pression 'SI' ('EGA' MOI 'PN'); 'FINSI'; *** captage de temperature 'SI' ('EGA' MOI 'TN'); 'FINSI'; 'FIN' BOU2; 'FIN' BOU1; 'FINPROC' TABC; ************************************************************************ ***************** FIN DEFINTION DES PROCEDURES ************************* ************************************************************************ ***************** REMPLISSAGE DE LA TABLE PGAZ ************************* * ***** Gas properties * PGAZ = 'TABLE' ; * Polynomial degree of specific heats PGAZ . 'NORD' = 4 ; * Species explicitly treated in the Euler Equations * Species non explicitly treated PGAZ . 'ESPNEULE' = 'N2 '; * Single gases properties PGAZ . 'H2 ' = 'TABLE' ; PGAZ . 'H2O ' = 'TABLE' ; PGAZ . 'N2 ' = 'TABLE' ; PGAZ . 'O2 ' = 'TABLE' ; * R (J/Kg/K) mH2 = 2. '*' 1.00797E-3 ; mo2 = 2. '*' 15.9994E-3 ; mH2O = mh2 '+' (0.5 '*' mo2) ; mN2 = 2 '*' 14.0067E-3 ; RGAS = 8.31441 ; PGAZ . 'H2 ' . 'R' = RGAS '/' mh2 ; PGAZ . 'H2O ' . 'R' = RGAS '/' mh2o ; PGAZ . 'N2 ' . 'R' = RGAS '/' mn2 ; PGAZ . 'O2 ' . 'R' = RGAS '/' mo2 ; * Polynomials regressions coefficients -2.37281455E-07 1.84701105E-11 ; -1.82753232E-08 2.44485692E-12 ; 8.78233606E-09 -3.05514485E-13 ; 2.33636971E-08 -1.53304905E-12; * Formation enthalpies at 0K (J/Kg)) * * h_i(0K) = h_i(T0) '-' \int_0^{T0} cp_i(x) dx * = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0) * PGAZ . 'H2 ' . 'H0K' = -4.195D6 ; PGAZ . 'H2O ' . 'H0K' = -1.395D7 ; PGAZ . 'N2 ' . 'H0K' = -2.953D5 ; PGAZ . 'O2 ' . 'H0K' = -2.634D5 ; *'OPTION' 'SAUVER' 'FORMAT' 'pgaz.sauv' ; *'SAUVER' 'FORMAT' PGAZ ; *'FIN'; *** ------------------------------------------------- *** *** Fichier correspondant a la definition du maillage *** *** ------------------------------------------------- *** ************ * MAILLAGE * ************ raf = 1 ; * * * Maillage du canal RUT sans obstacles * * *NX1 = 10 ; NX2 = 20 ; NX3 = 106 ; NX4 = 180 ; *DINI1 = 0.1 ; DFIN1 = 1. ; *NY1 = 15 ; NY2 = 25 ; NY3 = 10 ; NY4 = 13 ; X0 = 31.6 ; X1 = X0 ; Y0 = 0.0 ; Y1 = Y0 ; IGNI1 = 1. ; LC1 = 3. ; LK1 = 10.6 ; LC2 = 18. ; LC3 = 54. ; HK1 = 4.; HC = 2.3 ; *YJH = Y1 + 1. - ((HC - 1.)/(2. * NY4 - 1.));repositionnement *YJB = 1.5 - HK1 + ((1.5 - Y1) / ( 2. * NY1 - 1.)) ;repositionnement YJH = Y1 + 1. ; YJB = 1.5 - HK1 ; *---------------------------------------> Points du contour P1 = X1 Y1 ; P2 = (X1 + LC1) Y1 ; P3 = (X1 + LC1) (Y1 - HK1) ; P4 = (X1 + LC1 + LK1) (Y1 - HK1) ; P5 = (X1 + LC1 + LK1) Y1 ; PF1 = (X1 + LC1 + LK1 + LC2) Y1 ; PF2 = (X1 + LC1 + LK1 + LC2) (Y1 + HC) ; P6 = (X1 + LC1 + LK1 + LC2 + LC3) Y1 ; P7 = (X1 + LC1 + LK1 + LC2 + LC3) (Y1 + HC) ; P8 = (X1 + LC1 + LK1) (Y1 + HC) ; P9 = (X1 + LC1) (Y1 + HC) ; P10 = X1 (Y1 + HC) ; *---------------------------------------> Points pour les jauges J1 = X1 YJH ; J2 = (X1 + IGNI1) YJH ; J2B = (X1 + IGNI1) Y1 ; J2H = (X1 + IGNI1) (Y1 + HC) ; J3 = (X1 + LC1) YJH ; J4 = (X1 + LC1 + LK1) YJH ; J5 = (X1 + LC1 + LK1 + LC2) YJH ; J6 = (X1 + LC1 + LK1 + LC2 + LC3) YJH ; J7 = (X1 + LC1) YJB ; J8 = (X1 + LC1 + LK1) YJB ; * *- Lignes pour les sous-domaines * * MAILLAGE FIN *NX1 = 10 ; NX2 = 20 ; NX3 = 106 ; NX4 = 180 ; *DINI1 = 0.1 ; DFIN1 = 1. ; *NY1 = 15 ; NY2 = 25 ; NY3 = 10 ; NY4 = 13 ; * MAILLAGE TRES FIN *NX1 = 15 ; NX2 = 30 ; NX3 = 150 ; NX4 = 270 ; *DINI1 = 0.1 ; DFIN1 = 1. ; *NY1 = 22 ; NY2 = 37 ; NY3 = 15 ; NY4 = 19 ; * MAILLAGE TRES GROSSIER NX1 = 2 ; NX2 = 4 ; NX3 = 25 ; NX4 = 36 ; DINI1 = 0.5 ; DFIN1 = 5. ; NY1 = 3 ; NY2 = 5 ; NY3 = 2 ; NY4 = 3 ; * MAILLAGE GROSSIER *NX1 = 4 ; NX2 = 8 ; NX3 = 50 ; NX4 = 72 ; *DINI1 = 0.25 ; DFIN1 = 2.5 ; *NY1 = 6 ; NY2 = 10 ; NY3 = 4 ; NY4 = 6 ; *---------------------------------------> Contour externe *----------------------------------------> Jauges *----------------------------------------> Autres lignes internes * *- Pavage de chacune de zones identifiées * *---------------------------------------> Premier canal DOM1 = DC1 ET DC2 ; DOM2 = DC3 ET DC4 ; DOMC1 = DOM1 ET DOM2 ; *---------------------------------------> Deuxieme canal DOM3 = DC5 ET DC6 ; DOM4 = DC7 ET DC8 ; DOMC2 = DOM3 ET DOM4 ; *---------------------------------------> Canyon DOMK1 = DK1 ET DK2 ET DK3 ET DK4 ; *---------------------------------------> Maillage global DOMTOT = DOMC1 ET DOMC2 ET DOMK1 ; DDOM1 = DOM1 ; DDOM2 = DOM2 ET DOMC2 ET DOMK1 ; * *** PCEN1 Pour determiner PRES(PCEN,t) * XINIT0 = 0.33 ; *repositionnement *XINIT0 = 1.41 ; YINIT0 = 0.79 ; PCEL2 = (XINIT0 + 34.78) YJB ; PCEL4 = (XINIT0 + 41.29) YJB ; PCEL5 = (XINIT0 + 42.97) YJB ; PCEL6 = (XINIT0 + 44.84) YJB ; PCEL7 = (XINIT0 + 34.54) YJH ; PCEL8 = (XINIT0 + 36.79) YJH ; PCEL9 = (XINIT0 + 39.04) YJH ; PCEL10= (XINIT0 + 41.74) YJH ; PCEL11= (XINIT0 + 45.04) YJH ; PCEL12= (XINIT0 + 50.02) YJH ; PCEL13= (XINIT0 + 55.04) YJH ; PCEL14= (XINIT0 + 64.00) YJH ; * *PCEL2 = 35.11 YJB ; *PCEL4 = 41.62 YJB ; *PCEL5 = 43.30 YJB ; *PCEL6 = 45.17 YJB ; *PCEL7 = 34.87 YJH ; *PCEL8 = 37.12 YJH ; *PCEL9 = 39.37 YJH ; *PCEL10= 42.07 YJH ; *PCEL11= 45.37 YJH ; *PCEL12= 50.35 YJH ; *PCEL13= 55.37 YJH ; *PCEL14= 64.33 YJH ; LPOINTS = PCEL2 'ET' PCEL4 'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9 'ET' PCEL10 'ET' PCEL11 'ET' PCEL12 'ET' PCEL13 'ET' PCEL14 ; *'OPTI' 'SAUV' 'FORMAT' 'maillage_rut_fin.save'; *'OPTI' 'SAUV' 'FORMAT' 'maillage_rut_gros.save'; *'SAUV' 'FORMAT' ; *'FIN'; ******* FICHIER MAITRE ************************************* *** RECUPERATION DES OBJETS MODELES + POINTS DES CAPTEURS *** *'OPTI' 'REST' 'FORMAT' 'maillage_rut_gros.save'; *'OPTI' 'REST' 'FORMAT' 'maillage_rut_fin.save'; *'REST' 'FORMAT' ; * TABC = 'TABLE' 'CAPTEUR' ; TABC . 'LPOINTS' = PCEL2 'ET' PCEL4 'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9 'ET' PCEL10 'ET' PCEL11 'ET' PCEL12 'ET' PCEL13 'ET' PCEL14 ; 'C8' 'C9' 'C10' 'C11' 'C12' 'C13' 'C14' ; TABC . 'FREQ' = 1 ; *** DEFINITION DES GAZ EN PRESENCE + LEURS PROPRIETES ******** *********************** **** LA TABLE PGAZ **** *********************** *'OPTI' 'REST' 'FORMAT' 'pgaz.sauv'; *'REST' 'FORMAT' ; * *** Fin PGAZ * *** DEFINITION DES CONDITIONS INITIALES ************************* CONDINI = 'TABLE' ; CONDINI . 'NBDOM' = 2 ; * 1.1) DOM1 = gauche * 1.2) DOM2 = droit CONDINI . 'DOM1' = 'TABLE' ; CONDINI . 'DOM2' = 'TABLE' ; * ZONE GAUCHE tg = 3000.0 ; pg =20.00D5 ; xh2sec = 0.250 ; xh2o = 0.250 ; xh2 = xh2sec * (1.D0 - xh2o) ; xo2 = 0.21 * (1.D0 - xh2 - xh2o) ; xn2 = 1.0D0 - (xh2 + xo2 + xh2O) ; uxg = 0.0 ; uyg = 0.0 ; CONDINI . 'DOM1' .'MAIL' = DDOM1; CONDINI . 'DOM1' .'T' = tg ; CONDINI . 'DOM1' .'P' = pg ; CONDINI . 'DOM1' .'V' = (uxg uyg) ; CONDINI . 'DOM1' .'XH2' = xh2 ; CONDINI . 'DOM1' .'XO2' = xo2 ; CONDINI . 'DOM1' .'XN2' = xn2 ; CONDINI . 'DOM1' .'XH2O' = xh2o ; * ZONE DROITE td = 363. ; und = 0.0D0 ; utd = 0.0D0 ; pd = 1.0D5 ; CONDINI . 'DOM2' .'MAIL' = DDOM2; CONDINI . 'DOM2' .'T' = td ; CONDINI . 'DOM2' .'P' = pd ; CONDINI . 'DOM2' .'V' = (und utd) ; CONDINI . 'DOM2' .'XH2' = xh2 ; CONDINI . 'DOM2' .'XO2' = xo2 ; CONDINI . 'DOM2' .'XN2' = xn2 ; CONDINI . 'DOM2' .'XH2O' = xh2o ; *** Methodes possibles : * * 'VLH ' CONDINI . 'METO' = 'VLH ' ; CONDINI . 'MAILTOT' = DOMTOT ; CONDINI . 'TINI' = 0.0 ; CONDINI . 'TFINA' = 8.D-3 ; CONDINI . 'ORDREESP' = 2 ; * Maillage fin **CONDINI . 'DTIMP' = 1.D-5 ; * Maillage tres grossier CONDINI . 'DTIMP' = 3.D-5 ; *CONDINI . 'GRAPH' = VRAI ; CONDINI . 'GRAPH' = FAUX ; 'OPTI' 'TRACE' 'PSC' ; *** FIN remplissage de la table definissant les conditions initiales *** preparation des equations et initialisation des champoints *** recuperation de la table de calculs RV = PREPADET CONDINI PGAZ TABC ; MDOMT = RV.'GEO'.'MQUATOT' ; DOMT1 = RV.'GEO'.'MAILTOT'; X1 = 'COORD' 1 DOMT1; TDOMT1 . 'PRECONDI' = 1 ; * Table perso pour le trace de courbes au cours du temps *RV.'PRCPERSO' = 'TABLE' 'PERSO' ; *RV.'PRCPERSO'.'NOMPROC' = 'MOT' PERSO ; *RV.'PRCPERSO'.'ARG1' = $DOMT1 ; *RV.'PRCPERSO'.'ARG2' = 1 ; * Maillage fin **RV.'PRCPERSO'.'ARG3' = 15 ; * Maillage grossier *RV.'PRCPERSO'.'ARG3' = 5 ; *** execution du calcul EXEX RV ; *** relancement d un calcul *RV . 'PASDETPS' . 'TFINA' = 6.D-3 ; * EXEX RV ; * Maillage fin **'OPTION' 'SAUV' 'main_rut_fin.sauv' ; * Maillage grossier *'OPTION' 'SAUV' 'main_rut_gros.sauv' ; *'SAUV' ; *************** TESTS de non-regression ******************************** TPS = RV.'TABC'.'LTPS'; CA = 'TABLE' ; CA . 1 = RV.'TABC'.'C7'.'PN'; CA . 2 = RV.'TABC'.'C8'.'PN'; CA . 3 = RV.'TABC'.'C9'.'PN'; CA . 4 = RV.'TABC'.'C10'.'PN'; TOL1 = 1.E-6 ; NN1 = 1 ; 'REPETER' BLO1 N1 ; 'SI' (T1 'EGA' T2 TOL1); NN1 = NN1 + 1 ; 'FINSI'; 'QUITTER' BLO1 ; 'FINSI'; 'FIN' BLO1 ; 'MESS' 'Pression de reference' ; 'LIST' PREF; 'MESS' 'Pression calculee' ; 'LIST' PCAL; * AB: Je definis PREFS pour calculer une erreur relative. PREFS = 'MINIMUM' PREF ; EPS_P = 1.E-5 ; ERROR = 0 ; **** TEST sur la pression 'SI' ((DIFF1 '/' PREFS) '>' EPS_P) ; 'MESS' 'Probleme sur la pression ' ; ERROR = ERROR + 1 ; 'FINSI'; **** Calcul de la vitesse de flamme X1 Y1 = COORD (RV.'TABC'.'RPOINTS'); P7 = POINT P0 5; P8 = POINT P0 6; P9 = POINT P0 7; P10 = POINT P0 8; DX8 = X8 - X7 ; DX9 = X9 - X8 ; DX10 = X10 - X9 ; DT1 = T2 - T1 ; US1 = DX1 / DT1 ; 'SI' (ERROR '>' 0) ; 'FINSI'; 'FIN';
© Cast3M 2003 - Tous droits réservés.
Mentions légales