1 : C FLACR3 SOURCE BECC 05/06/09 21:15:07 5106 2 : SUBROUTINE FLACR3(EPSCSI,EPSI,DELTAT,MCEN,MFACEL,IRC,IYC,IYINIT, 3 : & IYFINA,IVCAR,IDX,MLRMAS,MLRH0K,MLRECO,ICHRET,ICHRY) 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : FLACR3 9 : C 10 : C DESCRIPTION : CREBCOM: modele non-homogene 11 : C voir FLACR2 12 : C 13 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) 14 : C 15 : C AUTEUR : A. BECCANTINI, DM2S/SFME/LTMF 16 : C 17 : C************************************************************************ 18 : C 19 : C INPUTS 20 : C 21 : C EPSCSI : (REAL*8), parametre pour controler si CSI=DY/(YF - YI) 22 : C \in(0,1). Si CSI1+EPSCSI un message de 23 : C warnin est donné. 24 : C 25 : C EPSI : parametre du critère CREBCOM (REAL*8) 26 : C 27 : C DELTAT : pas de temps (REAL*8) 28 : C 29 : C MCEN : pointeur du MELEME contenant les centres des ELTs 30 : C 31 : C MFACEL : pointeur du MELEME contenant la correspondence 32 : C CENTRE-FACE-CENTRE 33 : C 34 : C IRC : pointeur du CHPOINT contenant la masse volumique. 35 : C 36 : C IYC : pointeur du CHPOINT contenant les fractions massiques 37 : C 38 : C IYINIT : pointeur du CHPOINT contenant la fraction massique initiale 39 : C de la premiere composante de IYC; 40 : C 41 : C IYFINA : pointeur du CHPOINT contenant la fraction massique finale de 42 : C la premiere composante de IYC ; 43 : C 44 : C IVCAR : pointeur du CHPOINT contenant la vitesse caractéristique de 45 : C la flamme 46 : C 47 : C IDX : pointeur du CHPOINT contenant la diemnsion de la maille 48 : C 49 : C MLRMAS : pointeur du LISTREEL contenant les masses molaires 50 : C 51 : C MLRH0K : pointeur du LISTREEL contenant les energies des formation à 52 : C 0K 53 : C 54 : C MLRCOE : pointeur du LISTREEL contenant les coeff. stoch. de la 55 : C reaction chimique 56 : C 57 : C OUTPUTS : 58 : C 59 : C ICHRET : pointeur de l'increment de l'energie totale par unité de 60 : C volume 61 : C 62 : C ICHRY : pointeur de l'increment des densités massiques 63 : C 64 : C************************************************************************ 65 : C 66 : C HISTORIQUE (Anomalies et modifications éventuelles) 67 : C 68 : C HISTORIQUE : 69 : C 70 : C 71 : C************************************************************************ 72 : C 73 : IMPLICIT INTEGER(I-N) 74 : -INC CCOPTIO 75 : -INC SMLREEL 76 : POINTEUR MLRECO.MLREEL, MLRMAS.MLREEL, MLRH0K.MLREEL 77 : -INC SMELEME 78 : POINTEUR MELCEN.MELEME, MELEFL.MELEME 79 : -INC SMCHPOI 80 : INTEGER N, NC 81 : POINTEUR MPOVRO.MPOVAL, MPOVAY.MPOVAL,MPOVYI.MPOVAL, 82 : & MPOVYF.MPOVAL, MPOVC.MPOVAL, MPODRE.MPOVAL, MPODRY.MPOVAL, 83 : & MPOCSI.MPOVAL,MPOCRI.MPOVAL,MPOVDX.MPOVAL 84 : -INC SMLENTI 85 : POINTEUR MLECEN.MLENTI 86 : C 87 : C 88 : C**** Variables de COOPTIO 89 : C 90 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 91 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 92 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 93 : C & ,IECHO, IIMPI, IOSPI 94 : C & ,IDIM 95 : C & ,MCOORD 96 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 97 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 98 : C & ,NORINC,NORVAL,NORIND,NORVAD 99 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI 100 : C 101 : C**** Les variables 102 : C 103 : INTEGER MCEN,MFACEL,IRC,IYC,IYINIT,IYFINA,IVCAR,ICHRET,ICHRY 104 : $ ,IGEOM,ICEN,NCEN,NFAC,NLCF,NGCEG,NGCED,NLCEG,NLCED 105 : & ,IESP,NESP,IDX 106 : REAL*8 EPSI, DELTAT, Y1F, Y1I, Y1, VCSIG, VCSID, VCSI2G, VCSI2D 107 : & ,EPS12, DY, DYMAX, DCSI, DCSI1, DY1, RHO, DY2, YMAX, DYMAX1 108 : & ,ERRTOL, EPSCSI, CSIG 109 : CHARACTER*8 TYPE 110 : PARAMETER(ERRTOL=1.0D-6) 111 : C 112 : NESP=MLRH0K.PROG(/1) 113 : MELCEN=MCEN 114 : MELEFL=MFACEL 115 : C 116 : C**** KRIPAD pour la correspondance global/local de centre 117 : C 117 : CALL KRIPAD(MELCEN,MLECEN) 118 : IF(IERR .NE. 0)GOTO 9999 119 : C SEGINI MLECEN 120 : SEGACT MELCEN 121 : NCEN=MELCEN.NUM(/2) 122 : SEGACT MELEFL 123 : NFAC=MELEFL.NUM(/2) 124 : C 125 : C**** Lectures de CHPOINT 126 : C 126 : CALL LICHT(IRC,MPOVRO,TYPE,IGEOM) 127 : C SEGACT MPOVRO*MOD 127 : CALL LICHT(IYC,MPOVAY,TYPE,IGEOM) 128 : C SEGACT MPOVAY*MOD 128 : CALL LICHT(IYINIT,MPOVYI,TYPE,IGEOM) 129 : C SEGACT MPOVYI*MOD 129 : CALL LICHT(IYFINA,MPOVYF,TYPE,IGEOM) 130 : C SEGACT MPOVYF*MOD 130 : CALL LICHT(IVCAR,MPOVC,TYPE,IGEOM) 131 : C SEGACT MPOVYF*MOD 131 : CALL LICHT(IDX,MPOVDX,TYPE,IGEOM) 132 : C SEGACT MPOVDX*MOD 132 : CALL LICHT(ICHRET,MPODRE,TYPE,IGEOM) 133 : C SEGACT MPODRE*MOD 133 : CALL LICHT(ICHRY,MPODRY,TYPE,IGEOM) 134 : C SEGACT MPODRY*MOD 135 : C 136 : C**** Creation du MPOVAL qui contient csi 137 : C Creation du MPOVAL du critere 138 : C 139 : N=NCEN 140 : NC=1 141 : SEGINI MPOCSI 142 : SEGINI MPOCRI 143 : C 144 : C**** Calcul de DYMAX pour la premiere espece 145 : C Controle de la densité positive 146 : C Controle des fractions massiques (Yi>0, sum_i Y_i < 1) 147 : C 148 : DYMAX = 0.0D0 149 : YMAX=0.0D0 150 : DO ICEN=1,NCEN,1 151 : Y1F = MPOVYF.VPOCHA(ICEN,1) 152 : Y1I = MPOVYI.VPOCHA(ICEN,1) 153 : DY = Y1F - Y1I 154 : DYMAX=MAX(ABS(DY),DYMAX) 155 : Y1=MPOVAY.VPOCHA(ICEN,1) 156 : IF(Y1 .LT. 0.0D0)THEN 157 : WRITE(IOIMP,*) 'CHPO2 < 0 ???' 157 : CALL ERREUR(21) 158 : GOTO 9999 159 : ELSE 160 : YMAX=Y1 161 : ENDIF 162 : RHO=MPOVRO.VPOCHA(ICEN,1) 163 : IF(RHO .LT. 0)THEN 164 : WRITE(IOIMP,*) 'CHPO1 < 0 ???' 164 : CALL ERREUR(21) 165 : GOTO 9999 166 : ENDIF 167 : RHO=MPOVYI.VPOCHA(ICEN,1) 168 : IF((RHO .LT. 0) .OR. (RHO .GT. 1))THEN 169 : WRITE(IOIMP,*) 'CHPO3 < 0 ou CHPO3 > 1???' 169 : CALL ERREUR(21) 170 : GOTO 9999 171 : ENDIF 172 : RHO=MPOVYF.VPOCHA(ICEN,1) 173 : IF((RHO .LT. 0.0D0) .OR. (RHO .GT. 1.0D0))THEN 174 : WRITE(IOIMP,*) 'CHPO4 < 0 ou CHPO4 > 1???' 174 : CALL ERREUR(21) 175 : GOTO 9999 176 : ENDIF 177 : RHO=MPOVC.VPOCHA(ICEN,1) 178 : IF((RHO .LT. 0))THEN 179 : WRITE(IOIMP,*) 'CHPO5 < 0' 179 : CALL ERREUR(21) 180 : GOTO 9999 181 : ENDIF 182 : DO IESP=2,NESP,1 183 : Y1=MPOVAY.VPOCHA(ICEN,IESP) 184 : IF(Y1 .LT. 0.0D0)THEN 185 : WRITE(IOIMP,*) 'CHPO2 < 0 ???' 185 : CALL ERREUR(21) 186 : GOTO 9999 187 : ELSE 188 : YMAX=YMAX+Y1 189 : ENDIF 190 : ENDDO 191 : IF(YMAX .GT. 1.0D0)THEN 192 : WRITE(IOIMP,*) 'sum CHPO2 > 1 ???' 192 : CALL ERREUR(21) 193 : GOTO 9999 194 : ENDIF 195 : ENDDO 196 : C 197 : C**** Calcul de CSI \in (0,1) 198 : C 199 : DO ICEN=1,NCEN,1 200 : Y1F = MPOVYF.VPOCHA(ICEN,1) 201 : Y1I = MPOVYI.VPOCHA(ICEN,1) 202 : DY = Y1F - Y1I 203 : IF(ABS(DY) .LE. (ERRTOL*DYMAX))THEN 204 : C 205 : C********** Pas de combustion en cette region 206 : C 207 : MPOCSI.VPOCHA(ICEN,1)=0.0D0 208 : C 209 : ELSE 210 : Y1 = MPOVAY.VPOCHA(ICEN,1) 211 : VCSIG = (Y1 - Y1I) / DY 212 : C 213 : C********** On n'accepte pas csi > 1.1 or csi < -0.1 214 : C 215 : IF((VCSIG .GT. (1.0D0+EPSCSI)) .OR. 216 : & (VCSIG .LT. (-1*EPSCSI)))THEN 217 : WRITE(IOIMP,*) 'Progress variable = ???' 218 : C 21 2 219 : C Données incompatibles 219 : CALL ERREUR(21) 220 : GOTO 9999 221 : ELSEIF(VCSIG .GT. 1.0D0)THEN 222 : VCSIG = 1.0D0 223 : ELSEIF(VCSIG .LT. 0.0D0)THEN 224 : VCSIG = 0.0D0 225 : ENDIF 226 : MPOCSI.VPOCHA(ICEN,1)= VCSIG 227 : ENDIF 228 : ENDDO 229 : C 230 : C**** Le critere CREBCOM 231 : C 232 : DO NLCF = 1, NFAC, 1 233 : C 234 : C******* NLCF = numero local du centre de facel 235 : C NGCEG = numero global du centre ELT "gauche" 236 : C NLCEG = numero local du centre ELT "gauche" 237 : C NGCED = numero global du centre ELT "droite" 238 : C NLCED = numero local du centre ELT "droite" 239 : C 240 : NGCEG = MELEFL.NUM(1,NLCF) 241 : NGCED = MELEFL.NUM(3,NLCF) 242 : NLCEG = MLECEN.LECT(NGCEG) 243 : NLCED = MLECEN.LECT(NGCED) 244 : C 245 : VCSIG=MPOCSI.VPOCHA(NLCEG,1) 246 : VCSID=MPOCSI.VPOCHA(NLCED,1) 247 : VCSI2G=VCSIG*VCSIG 248 : VCSI2D=VCSID*VCSID 249 : C 250 : IF(NLCEG .EQ. NLCED)THEN 251 : C 252 : C********** Murs 253 : C 254 : MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) + (0.5D0 * 255 : & VCSI2D) 256 : C 257 : ELSE 258 : C 259 : MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) + 260 : & (VCSI2D - (0.5D0 * VCSI2G)) 261 : MPOCRI.VPOCHA(NLCED,1)=MPOCRI.VPOCHA(NLCED,1) + 262 : & (VCSI2G - (0.5D0 * VCSI2D)) 263 : C 264 : ENDIF 265 : ENDDO 266 : C 267 : C**** Calcul des increments de l'energie et des densités 268 : C 269 : EPS12 = EPSI * EPSI 270 : DO ICEN = 1, NCEN, 1 271 : CSIG = MPOCSI.VPOCHA(ICEN,1) 272 : VCSIG = MPOCRI.VPOCHA(ICEN,1) 273 : C 274 : C******* In 2D, contribution of the ideal upper and lower cells 275 : C 276 : IF(IDIM .EQ. 2) VCSIG = VCSIG + (CSIG * CSIG) 277 : IF(VCSIG .GE. (EPS12*(1.0D0+ERRTOL)))THEN 278 : C 279 : C********** Il y a combustion 280 : C 281 : DCSI=MPOVC.VPOCHA(ICEN,1)*DELTAT/MPOVDX.VPOCHA(ICEN,1) 282 : DCSI1=1.0D0-CSIG 283 : DCSI=MIN(DCSI,DCSI1) 284 : DY1=(MPOVYF.VPOCHA(ICEN,1)-MPOVYI.VPOCHA(ICEN,1))*DCSI 285 : C 286 : C********** On force la positivité de MPOVAY.VPOCHA(ICEN,1) 287 : C 288 : IF(DY1 .GT. 0.0D0)THEN 289 : DY2 = 1.0D0 - MPOVAY.VPOCHA(ICEN,1) 290 : ELSE 291 : DY2 = MPOVAY.VPOCHA(ICEN,1) 292 : ENDIF 293 : DY1=SIGN(MIN(ABS(DY2),ABS(DY1)),DY1) 294 : C 295 : RHO=MPOVRO.VPOCHA(ICEN,1) 296 : MPODRY.VPOCHA(ICEN,1)=DY1*RHO 297 : MPODRE.VPOCHA(ICEN,1)=-1.0D0*DY1*RHO*MLRH0K.PROG(1) 298 : DY1=DY1/(MLRMAS.PROG(1)*MLRECO.PROG(1)) 299 : DO IESP=2,NESP,1 300 : DY=DY1*(MLRMAS.PROG(IESP)*MLRECO.PROG(IESP)) 301 : DYMAX1=ABS(DYMAX/(MLRMAS.PROG(1)*MLRECO.PROG(1))* 302 : & (MLRMAS.PROG(IESP)*MLRECO.PROG(IESP))) 303 : IF(DY .GT. 0)THEN 304 : DY2=1.0D0 - MPOVAY.VPOCHA(ICEN,IESP) 305 : ELSE 306 : DY2=MPOVAY.VPOCHA(ICEN,IESP) 307 : ENDIF 308 : C 309 : C************* On force la positivité de MPOVAY.VPOCHA(ICEN,IESP) 310 : C 311 : IF((ABS(DY)-ABS(DY2)).GE.(0.5D0*DYMAX1))THEN 312 : WRITE(IOIMP,*) 'CHPO2, CHPO3, CHPO4 = ???' 313 : C Données incompatibles 313 : CALL ERREUR(21) 314 : GOTO 9999 315 : ELSE 316 : DY=SIGN(MIN(ABS(DY2),ABS(DY)),DY) 317 : ENDIF 318 : MPODRY.VPOCHA(ICEN,IESP)=DY*RHO 319 : MPODRE.VPOCHA(ICEN,1)=MPODRE.VPOCHA(ICEN,1)-(DY*RHO 320 : $ *MLRH0K.PROG(IESP)) 321 : ENDDO 322 : ENDIF 323 : ENDDO 324 : C 325 : SEGDES MELCEN 326 : SEGDES MELEFL 327 : SEGSUP MLECEN 328 : C 329 : SEGDES MPOVRO 330 : SEGDES MPOVAY 331 : SEGDES MPOVYI 332 : SEGDES MPOVYF 333 : SEGDES MPOVC 334 : SEGDES MPOVDX 335 : SEGDES MPODRE 336 : SEGDES MPODRY 337 : SEGSUP MPOCSI 338 : SEGSUP MPOCRI 339 : C 340 : 9999 RETURN 341 : END 342 : 343 : 344 : 345 : 346 : 347 : 348 : 349 : 350 : 351 : 352 : 353 : 354 :
© Cast3M 2003 - All rights reserved.
Disclaimer