C FLACR3 SOURCE CB215821 20/11/25 13:29:00 10792 SUBROUTINE FLACR3(EPSCSI,EPSI,DELTAT,MCEN,MFACEL,IRC,IYC,IYINIT, & IYFINA,IVCAR,IDX,MLRMAS,MLRH0K,MLRECO,ICHRET,ICHRY) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : FLACR3 C C DESCRIPTION : CREBCOM: modele non-homogene C voir FLACR2 C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) C C AUTEUR : A. BECCANTINI, DM2S/SFME/LTMF C C************************************************************************ C C INPUTS C C EPSCSI : (REAL*8), parametre pour controler si CSI=DY/(YF - YI) C \in(0,1). Si CSI1+EPSCSI un message de C warnin est donné. C C EPSI : parametre du critère CREBCOM (REAL*8) C C DELTAT : pas de temps (REAL*8) C C MCEN : pointeur du MELEME contenant les centres des ELTs C C MFACEL : pointeur du MELEME contenant la correspondence C CENTRE-FACE-CENTRE C C IRC : pointeur du CHPOINT contenant la masse volumique. C C IYC : pointeur du CHPOINT contenant les fractions massiques C C IYINIT : pointeur du CHPOINT contenant la fraction massique initiale C de la premiere composante de IYC; C C IYFINA : pointeur du CHPOINT contenant la fraction massique finale de C la premiere composante de IYC ; C C IVCAR : pointeur du CHPOINT contenant la vitesse caractéristique de C la flamme C C IDX : pointeur du CHPOINT contenant la diemnsion de la maille C C MLRMAS : pointeur du LISTREEL contenant les masses molaires C C MLRH0K : pointeur du LISTREEL contenant les energies des formation à C 0K C C MLRCOE : pointeur du LISTREEL contenant les coeff. stoch. de la C reaction chimique C C OUTPUTS : C C ICHRET : pointeur de l'increment de l'energie totale par unité de C volume C C ICHRY : pointeur de l'increment des densités massiques C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : C C C************************************************************************ C IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC SMLREEL POINTEUR MLRECO.MLREEL, MLRMAS.MLREEL, MLRH0K.MLREEL -INC SMELEME POINTEUR MELCEN.MELEME, MELEFL.MELEME -INC SMCHPOI INTEGER N, NC POINTEUR MPOVRO.MPOVAL, MPOVAY.MPOVAL,MPOVYI.MPOVAL, & MPOVYF.MPOVAL, MPOVC.MPOVAL, MPODRE.MPOVAL, MPODRY.MPOVAL, & MPOCSI.MPOVAL,MPOCRI.MPOVAL,MPOVDX.MPOVAL -INC SMLENTI POINTEUR MLECEN.MLENTI C C C**** Variables de COOPTIO C C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES C & ,IECHO, IIMPI, IOSPI C & ,IDIM C & ,MCOORD C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU C & ,NORINC,NORVAL,NORIND,NORVAD C & ,NUCROU, IPSAUV, IFICLE, IPREFI C C**** Les variables C INTEGER MCEN,MFACEL,IRC,IYC,IYINIT,IYFINA,IVCAR,ICHRET,ICHRY $ ,IGEOM,ICEN,NCEN,NFAC,NLCF,NGCEG,NGCED,NLCEG,NLCED & ,IESP,NESP,IDX REAL*8 EPSI, DELTAT, Y1F, Y1I, Y1, VCSIG, VCSID, VCSI2G, VCSI2D & ,EPS12, DY, DYMAX, DCSI, DCSI1, DY1, RHO, DY2, YMAX, DYMAX1 & ,ERRTOL, EPSCSI, CSIG CHARACTER*8 TYPE PARAMETER(ERRTOL=1.0D-6) C NESP=MLRH0K.PROG(/1) MELCEN=MCEN MELEFL=MFACEL C C**** KRIPAD pour la correspondance global/local de centre C CALL KRIPAD(MELCEN,MLECEN) IF(IERR .NE. 0)GOTO 9999 C SEGINI MLECEN SEGACT MELCEN NCEN=MELCEN.NUM(/2) SEGACT MELEFL NFAC=MELEFL.NUM(/2) C C**** Lectures de CHPOINT C CALL LICHT(IRC,MPOVRO,TYPE,IGEOM) C SEGACT MPOVRO*MOD CALL LICHT(IYC,MPOVAY,TYPE,IGEOM) C SEGACT MPOVAY*MOD CALL LICHT(IYINIT,MPOVYI,TYPE,IGEOM) C SEGACT MPOVYI*MOD CALL LICHT(IYFINA,MPOVYF,TYPE,IGEOM) C SEGACT MPOVYF*MOD CALL LICHT(IVCAR,MPOVC,TYPE,IGEOM) C SEGACT MPOVYF*MOD CALL LICHT(IDX,MPOVDX,TYPE,IGEOM) C SEGACT MPOVDX*MOD CALL LICHT(ICHRET,MPODRE,TYPE,IGEOM) C SEGACT MPODRE*MOD CALL LICHT(ICHRY,MPODRY,TYPE,IGEOM) C SEGACT MPODRY*MOD C C**** Creation du MPOVAL qui contient csi C Creation du MPOVAL du critere C N=NCEN NC=1 SEGINI MPOCSI SEGINI MPOCRI C C**** Calcul de DYMAX pour la premiere espece C Controle de la densité positive C Controle des fractions massiques (Yi>0, sum_i Y_i < 1) C DYMAX = 0.0D0 YMAX=0.0D0 DO ICEN=1,NCEN,1 Y1F = MPOVYF.VPOCHA(ICEN,1) Y1I = MPOVYI.VPOCHA(ICEN,1) DY = Y1F - Y1I DYMAX=MAX(ABS(DY),DYMAX) Y1=MPOVAY.VPOCHA(ICEN,1) IF(Y1 .LT. 0.0D0)THEN WRITE(IOIMP,*) 'CHPO2 < 0 ???' CALL ERREUR(21) GOTO 9999 ELSE YMAX=Y1 ENDIF RHO=MPOVRO.VPOCHA(ICEN,1) IF(RHO .LT. 0)THEN WRITE(IOIMP,*) 'CHPO1 < 0 ???' CALL ERREUR(21) GOTO 9999 ENDIF RHO=MPOVYI.VPOCHA(ICEN,1) IF((RHO .LT. 0) .OR. (RHO .GT. 1))THEN WRITE(IOIMP,*) 'CHPO3 < 0 ou CHPO3 > 1???' CALL ERREUR(21) GOTO 9999 ENDIF RHO=MPOVYF.VPOCHA(ICEN,1) IF((RHO .LT. 0.0D0) .OR. (RHO .GT. 1.0D0))THEN WRITE(IOIMP,*) 'CHPO4 < 0 ou CHPO4 > 1???' CALL ERREUR(21) GOTO 9999 ENDIF RHO=MPOVC.VPOCHA(ICEN,1) IF((RHO .LT. 0))THEN WRITE(IOIMP,*) 'CHPO5 < 0' CALL ERREUR(21) GOTO 9999 ENDIF DO IESP=2,NESP,1 Y1=MPOVAY.VPOCHA(ICEN,IESP) IF(Y1 .LT. 0.0D0)THEN WRITE(IOIMP,*) 'CHPO2 < 0 ???' CALL ERREUR(21) GOTO 9999 ELSE YMAX=YMAX+Y1 ENDIF ENDDO IF(YMAX .GT. 1.0D0)THEN WRITE(IOIMP,*) 'sum CHPO2 > 1 ???' CALL ERREUR(21) GOTO 9999 ENDIF ENDDO C C**** Calcul de CSI \in (0,1) C DO ICEN=1,NCEN,1 Y1F = MPOVYF.VPOCHA(ICEN,1) Y1I = MPOVYI.VPOCHA(ICEN,1) DY = Y1F - Y1I IF(ABS(DY) .LE. (ERRTOL*DYMAX))THEN C C********** Pas de combustion en cette region C MPOCSI.VPOCHA(ICEN,1)=0.0D0 C ELSE Y1 = MPOVAY.VPOCHA(ICEN,1) VCSIG = (Y1 - Y1I) / DY C C********** On n'accepte pas csi > 1.1 or csi < -0.1 C IF((VCSIG .GT. (1.0D0+EPSCSI)) .OR. & (VCSIG .LT. (-1*EPSCSI)))THEN WRITE(IOIMP,*) 'Progress variable = ???' C 21 2 C Données incompatibles CALL ERREUR(21) GOTO 9999 ELSEIF(VCSIG .GT. 1.0D0)THEN VCSIG = 1.0D0 ELSEIF(VCSIG .LT. 0.0D0)THEN VCSIG = 0.0D0 ENDIF MPOCSI.VPOCHA(ICEN,1)= VCSIG ENDIF ENDDO C C**** Le critere CREBCOM C DO NLCF = 1, NFAC, 1 C C******* NLCF = numero local du centre de facel C NGCEG = numero global du centre ELT "gauche" C NLCEG = numero local du centre ELT "gauche" C NGCED = numero global du centre ELT "droite" C NLCED = numero local du centre ELT "droite" C NGCEG = MELEFL.NUM(1,NLCF) NGCED = MELEFL.NUM(3,NLCF) NLCEG = MLECEN.LECT(NGCEG) NLCED = MLECEN.LECT(NGCED) C VCSIG=MPOCSI.VPOCHA(NLCEG,1) VCSID=MPOCSI.VPOCHA(NLCED,1) VCSI2G=VCSIG*VCSIG VCSI2D=VCSID*VCSID C IF(NLCEG .EQ. NLCED)THEN C C********** Murs C MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) + (0.5D0 * & VCSI2D) C ELSE C MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) + & (VCSI2D - (0.5D0 * VCSI2G)) MPOCRI.VPOCHA(NLCED,1)=MPOCRI.VPOCHA(NLCED,1) + & (VCSI2G - (0.5D0 * VCSI2D)) C ENDIF ENDDO C C**** Calcul des increments de l'energie et des densités C EPS12 = EPSI * EPSI DO ICEN = 1, NCEN, 1 CSIG = MPOCSI.VPOCHA(ICEN,1) VCSIG = MPOCRI.VPOCHA(ICEN,1) C C******* In 2D, contribution of the ideal upper and lower cells C IF(IDIM .EQ. 2) VCSIG = VCSIG + (CSIG * CSIG) IF(VCSIG .GE. (EPS12*(1.0D0+ERRTOL)))THEN C C********** Il y a combustion C DCSI=MPOVC.VPOCHA(ICEN,1)*DELTAT/MPOVDX.VPOCHA(ICEN,1) DCSI1=1.0D0-CSIG DCSI=MIN(DCSI,DCSI1) DY1=(MPOVYF.VPOCHA(ICEN,1)-MPOVYI.VPOCHA(ICEN,1))*DCSI C C********** On force la positivité de MPOVAY.VPOCHA(ICEN,1) C IF(DY1 .GT. 0.0D0)THEN DY2 = 1.0D0 - MPOVAY.VPOCHA(ICEN,1) ELSE DY2 = MPOVAY.VPOCHA(ICEN,1) ENDIF DY1=SIGN(MIN(ABS(DY2),ABS(DY1)),DY1) C RHO=MPOVRO.VPOCHA(ICEN,1) MPODRY.VPOCHA(ICEN,1)=DY1*RHO MPODRE.VPOCHA(ICEN,1)=-1.0D0*DY1*RHO*MLRH0K.PROG(1) DY1=DY1/(MLRMAS.PROG(1)*MLRECO.PROG(1)) DO IESP=2,NESP,1 DY=DY1*(MLRMAS.PROG(IESP)*MLRECO.PROG(IESP)) DYMAX1=ABS(DYMAX/(MLRMAS.PROG(1)*MLRECO.PROG(1))* & (MLRMAS.PROG(IESP)*MLRECO.PROG(IESP))) IF(DY .GT. 0)THEN DY2=1.0D0 - MPOVAY.VPOCHA(ICEN,IESP) ELSE DY2=MPOVAY.VPOCHA(ICEN,IESP) ENDIF C C************* On force la positivité de MPOVAY.VPOCHA(ICEN,IESP) C IF((ABS(DY)-ABS(DY2)).GE.(0.5D0*DYMAX1))THEN WRITE(IOIMP,*) 'CHPO2, CHPO3, CHPO4 = ???' C Données incompatibles CALL ERREUR(21) GOTO 9999 ELSE DY=SIGN(MIN(ABS(DY2),ABS(DY)),DY) ENDIF MPODRY.VPOCHA(ICEN,IESP)=DY*RHO MPODRE.VPOCHA(ICEN,1)=MPODRE.VPOCHA(ICEN,1)-(DY*RHO $ *MLRH0K.PROG(IESP)) ENDDO ENDIF ENDDO C SEGDES MELCEN SEGDES MELEFL SEGSUP MLECEN C SEGDES MPOVRO SEGDES MPOVAY SEGDES MPOVYI SEGDES MPOVYF SEGDES MPOVC SEGDES MPOVDX SEGDES MPODRE SEGDES MPODRY SEGSUP MPOCSI SEGSUP MPOCRI C 9999 RETURN END