1 : C DETO2 SOURCE CHAT 06/06/01 21:16:08 5450 2 : SUBROUTINE DETO2(X1,X2,X3,X4,P1,T1, 3 : 1 RAD,TAD,PAD,FACT,VCJ, 4 : 2 RZ,TZ,PZ, 5 : 3 RAICC,TAICC,PAICC, 6 : 4 AN,AH2,AH2O,AN2,AO2,TINT,IINT,NBC,NBT,IRET) 7 : C----------------------------------------------------------------------- 8 : C Calcul pour un mélange donné des états CJ, ZND et AICC 9 : C----------------------------------------------------------------------- 10 : C 11 : C--------------------------- 12 : C Parametres Entree/Sortie : 13 : C--------------------------- 14 : C 15 : C E/ X1 : Nombre de moles de H2 16 : C E/ X2 : Nombre de moles de O2 17 : C E/ X3 : Nombre de moles de H2O 18 : C E/ X4 : Nombre de moles de N2 19 : C E/ P1 : Pression du mélange initiale (Pa) 20 : C E/ T1 : Température du mélange initiale (K) 21 : C /S RAD : Densité pour l'état CJ (kg/m3) 22 : C /S TAD : Température pour l'état CJ (K) 23 : C /S PAD : Pression pour l'état CJ (Pa) 24 : C /S FACT : Taux de combustion pour l'état CJ (-) 25 : C /S VCJ : Vitesse de CJ (m/s) 26 : C /S RZ : Densité pour l'état ZND (kg/m3) 27 : C /S TZ : Température pour l'état ZND (K) 28 : C /S PZ : Pression pour l'état ZND (Pa) 29 : C /S RAICC : Densité pour l'état AICC (kg/m3) 30 : C /S TAICC : Température pour l'état AICC (K) 31 : C /S PAICC : Pression pour l'état AICC (Pa) 32 : C E/ AN(NBT) : Masse molaire des especes (en gramme) 33 : C E/ AH2(NBT) : Coeff du devt en T des propriétés de h2 34 : C E/ AH2O(NBT) : Coeff du devt en T des propriétés de h2o 35 : C E/ AN2(NBT) : Coeff du devt en T des propriétés de n2 36 : C E/ AO2(NBT) : Coeff du devt en T des propriétés de o2 37 : C E/ TINT : Température seuil (K) 38 : C E/ IINT : Nombre de coeff nécessaire à T donnée 39 : C E/ NBC : Nombre de constituants du mélange 40 : C E/ NBT : Nombre de coeff total 41 : C /S IRET : Indicateur de succès (0=OK, 1=Problème) 42 : C 43 : C----------------------------------------------------------------------- 44 : C 45 : C--------------------------- 46 : C Principe de l'algorithme : 47 : C--------------------------- 48 : C 49 : C 1) On part du point représentatif de l'état initial 50 : C 2) On suppose la combustion complète 51 : C 3) On calcule le point CJ itérativement en diminuant le volume 52 : C spécifique à partir des conditions suivantes : le point est sur 53 : C l'adiabatique d'hugoniot des produits et la tangente au point 54 : C passe par l'état initial (deux newtons imbriqués) 55 : C 4) On calcule le taux de production/destruction de H2 56 : C 5) Si équilibre, le point obtenu est le point de CJ (déto stable) ; 57 : C Sinon, on diminue le taux de combustion et on repart en 3) 58 : C (newton le plus externe) 59 : C 60 : C 6) On calcule ensuite le point ZND qui est sur l'hugoniot des gaz 61 : C frais et sur la droite de Rayleigh reliant le point initial au 62 : C point CJ (quand on aime newton on ne compte pas). 63 : C 64 : C 7) On calcule alors la vitesse de CJ à partir de la pente de la droite 65 : C de Rayleigh VCJ = ((PCJ - PI) / (VI - VCJ))**0.5 * VI 66 : C 67 : C----------------------------------------------------------------------- 68 : C 69 : C Langage : FORTRAN77 70 : C 71 : C Auteurs : F.DABBENE et E.STUDER 02/99 72 : C 73 : C----------------------------------------------------------------------- 74 : IMPLICIT INTEGER(I-N) 75 : IMPLICIT REAL*8 (A-H,O-Z) 76 : LOGICAL INIT,DROITE,GAUCHE,AICC 77 : DIMENSION X(4),WDOT(4),XC(4) 78 : DIMENSION AN(*),AH2(*),AH2O(*),AN2(*),AO2(*) 79 : DIMENSION NIT(9),NITER(9),TCX1(9),TCY1(9),TCX2(9),TCY2(9) 80 : DO ITC = 1, 9 81 : NIT(ITC)=0 82 : NITER(ITC)=1 83 : ENDDO 84 : C 85 : C- Récupération des données 86 : C 87 : IRET = 1 88 : X(1) = X1 89 : X(2) = X2 90 : X(3) = X3 91 : X(4) = X4 92 : C 93 : C- Fraction molaire de chaque constituant 94 : C 95 : XTOT = 0.D0 96 : DO 1 K = 1,NBC 97 : XTOT = XTOT + X(K) 98 : 1 CONTINUE 99 : DO 2 K = 1,NBC 100 : X(K) = X(K) / XTOT 101 : 2 CONTINUE 102 : C 103 : C- Enthalpie spécifique, densité et volume spécifique de l'état initial 104 : C 104 : CALL DETOH(T1,X,AN,AH2,AH2O,AN2,AO2,TINT,H1,IINT,NBC,NBT) 104 : CALL DETORO(P1,T1,X,AN,R1,NBC) 105 : V1 = 1.D0/R1 106 : C 107 : C- On verifie que le mélange est susceptible de déroner. 108 : C- Dans le cas contraire, on renvoie la condition initiale 109 : C 110 : IF ((X(1).LE.0.04D0) .OR. (X(3).GE.0.6D0)) THEN 111 : RAD = R1 112 : TAD = T1 113 : PAD = P1 114 : FACT = 0.D0 115 : VCJ = 0.D0 116 : RZ = R1 117 : TZ = T1 118 : PZ = P1 119 : RAICC = R1 120 : TAICC = T1 121 : PAICC = P1 122 : IRET = 0 123 : RETURN 124 : ENDIF 125 : C 126 : C- Initialisation du taux de combustion (1=complète) 127 : C 128 : FACT = 1.D0 129 : AICC = .TRUE. 130 : C--------------------------------------------------------------------- 131 : C- Retour du NEWTON sur le taux de combustion : la variable pivot est 132 : C- le taux de combustion. 133 : C--------------------------------------------------------------------- 134 : 70 CONTINUE 135 : C 136 : C- Evaluation des fractions molaires pour un taux de combustion égal à 137 : C- FACT et prise en compte d'un éventuel déficit en oxygène du mélange 138 : C 139 : XC(1) = (X(1) - FACT*X(1)) * P1/T1 140 : XC(2) = (X(2) - 0.5D0*FACT*X(1)) * P1/T1 141 : XC(3) = (X(3) + FACT*X(1)) * P1/T1 142 : XC(4) = X(4) * P1/T1 143 : IF (XC(2).LT.0.D0) THEN 144 : XC(2) = 0.D0 145 : XC(1) = (X(1) - 2.D0*X(2)) * P1/T1 146 : XC(3) = (X(3) + 2.D0*X(2)) * P1/T1 147 : ENDIF 148 : XTOT = 0.D0 149 : DO 3 K = 1,NBC 150 : XTOT = XTOT + XC(K) 151 : 3 CONTINUE 152 : DO 4 K = 1,NBC 153 : XC(K) = XC(K) / XTOT 154 : 4 CONTINUE 155 : C 156 : C- Initialisation de la densité R0 à la densité initiale 157 : C 158 : R0 = R1 159 : C--------------------------------------------------------------------- 160 : C- Retour du NEWTON sur le test d'appartenance du point initial à 161 : C- la tangente à l'hugoniot passant par le point CJ candidat : 162 : C- la variable pivot est la densité R0. 163 : C--------------------------------------------------------------------- 164 : 60 CONTINUE 165 : R2 = R0 166 : INIT=.TRUE. 167 : DROITE = .FALSE. 168 : GAUCHE = .FALSE. 169 : C--------------------------------------------------------------------- 170 : C- Retour afin de calculer les états à gauche et à droite du point 171 : C- (PC,TC,R2) susceptible d'etre l'etat de CJ afin d'évaluer la 172 : C- tangente à l'hugoniot. 173 : C--------------------------------------------------------------------- 174 : 50 CONTINUE 175 : PC = P1 176 : C--------------------------------------------------------------------- 177 : C- Retour du NEWTON permettant d'évaluer l'état d'équilibre (PC,TC) 178 : C- la densité étant fixée (=R2=R0) : on cherche (PC,TC) appartenant 179 : C- à la courbe d'hugoniot pour un taux de réaction égal à FACT. La 180 : C- variable pivot est la pression PC (de facon indireste TC). 181 : C--------------------------------------------------------------------- 182 : 20 CONTINUE 183 : WTM2 = XC(1)*AN(1) + XC(2)*AN(2) + XC(3)*AN(3) + XC(4)*AN(4) 184 : WTM2 = WTM2 * 1.D-3 185 : TC = WTM2*PC/(R2*8.314D0) 185 : CALL DETOH(TC,XC,AN,AH2,AH2O,AN2,AO2,TINT,H2,IINT,NBC,NBT) 186 : FX = 2.D0*(H2-H1)-((PC-P1)*((1.D0/R1)+(1.D0/R2))) 187 : X1 = 1.D0 188 : X2 = 5.D-5 188 : CALL NEWTON(PC,FX,X1,X2,1000,3,NFIN3,NIT,NITER,TCX1,TCY1,TCX2, 189 : $ TCY2) 190 : IF (NFIN3 .EQ. 0) GOTO 20 191 : IF (NFIN3 .EQ. 2) RETURN 192 : C 193 : C- Le point (PC,TC,R2) trouvé est l'état CJ à condition que la pente 194 : C- à l'hugoniot passe par le point initial (P1,T1,R1). Pour évaluer la 195 : C- pente, on calcule un état gauche et un état droit de (PC,TC,R2) en 196 : C- modifiant la densité de 1% (par défaut et par excès) 197 : C 198 : IF (GAUCHE) THEN 199 : PG = PC 200 : VG = 1.D0/R2 201 : GAUCHE = .FALSE. 202 : ENDIF 203 : IF (DROITE) THEN 204 : PD = PC 205 : VD = 1.D0/R2 206 : DROITE = .FALSE. 207 : GAUCHE = .TRUE. 208 : ENDIF 209 : IF (INIT) THEN 210 : PAD = PC 211 : TAD = TC 212 : VAD = 1.D0/R2 213 : C 214 : C- La première convergence donne l'état AICC puisque le point obtenu 215 : C- est le point sur l'hugoniot correspondant à un taux de combustion 216 : C- complète et à la densité initiale. 217 : C 218 : IF (AICC) THEN 219 : PAICC = PC 220 : TAICC = TC 221 : VAICC = 1.D0/R2 222 : AICC = .FALSE. 223 : ENDIF 224 : INIT = .FALSE. 225 : DROITE = .TRUE. 226 : ENDIF 227 : IF (DROITE) THEN 228 : R2 = R0*0.99D0 229 : GOTO 50 230 : ENDIF 231 : IF (GAUCHE) THEN 232 : R2 = R0*1.01D0 233 : GOTO 50 234 : ENDIF 235 : C 236 : C- Equation de la droite passant par les points (PTEST,V1) et (PC,VC) 237 : C- P = ALP*V + BCS 238 : C- On vérifie l'ecart entre PTEST et P1 ; en cas de non convergence 239 : C- on définit un nouveau R0 et on recommence la détermination d'un 240 : C- état CJ potentiel. 241 : C 242 : ALP = (PD-PG)/(VD-VG) 243 : BCS = PAD - (ALP*VAD) 244 : PTEST = ALP*V1+BCS 245 : FY = PTEST - P1 246 : X11 = 1.D0 247 : X21 = 5.D-5 247 : CALL NEWTON(R0,FY,X11,X21,1000,2,NFIN2,NIT,NITER,TCX1,TCY1, 248 : $ tcx2,tcy2) 249 : IF (NFIN2 .EQ. 0) GOTO 60 250 : IF (NFIN2 .EQ. 2) RETURN 251 : C 252 : C- XC, PAD, TAD et VAD caractérisent un état CJ qui est le bon à 253 : C- condition que l'équilibre chimique soit réalisé : On calcule 254 : C- le taux de production de H2 (i.e. de chaque espece puisque ici 255 : C- cinétique à une réaction) qui doit etre nul ; sinon, on diminue 256 : C- le taux de combustion. 257 : C 257 : CALL DETOWP(PAD,TAD,XC,AH2,AH2O,AO2,TINT,WDOT,IINT,NBC,NBT) 258 : FZ = WDOT(1) 259 : X13 =-0.001D0 260 : X23 = 1.D-20 260 : CALL NEWTON(FACT,FZ,X13,X23,1000,1,NFIN1,NIT,NITER,TCX1,TCY1, 261 : $ TCX2,TCY2) 262 : IF (NFIN1 .EQ. 0) GOTO 70 263 : IF (NFIN1 .EQ. 2) RETURN 264 : C 265 : C- Initialisation du volume spécifique de l'état ZND et de la pente 266 : C- de la droite de Rayleigh 267 : C 268 : VZ = V1/5.D0 269 : PENTE = (PAD-P1) / (VAD-V1) 270 : C--------------------------------------------------------------------- 271 : C- Retour du NEWTON du calcul ZND : on cherche le point d'intersection 272 : C- entre l'hugoniot passant par le point initial et la droite passant 273 : C- part CJ et le point initial. La variable pivot est le volume 274 : C- spécifique VZ. 275 : C--------------------------------------------------------------------- 276 : 80 CONTINUE 277 : PZ = P1 + (PENTE*(VZ-V1)) 278 : WTM1 = X(1)*AN(1) + X(2)*AN(2) + X(3)*AN(3) + X(4)*AN(4) 279 : WTM1 = WTM1 * 1.D-3 280 : TZ = WTM1*PZ*VZ / 8.314D0 280 : CALL DETOH(TZ,X,AN,AH2,AH2O,AN2,AO2,TINT,HZ,IINT,NBC,NBT) 281 : FQ =-2.D0 * (HZ-H1) + (PZ-P1)*(V1+VZ) 282 : X1 =-0.01D0 283 : X2 = 5.D-5 284 : X2 = 5.D-10 284 : CALL NEWTON(VZ,FQ,X1,X2,1000,1,NFIN,NIT,NITER,TCX1,TCY1, 285 : $ TCX2,TCY2) 286 : IF (NFIN .EQ. 0) GOTO 80 287 : IF (NFIN .EQ. 2) RETURN 288 : C 289 : C- Calcul des densités à partir des volumes spécifiques 290 : C 291 : RAD = 1.D0/VAD 292 : RAICC = 1.D0/VAICC 293 : RZ = 1.D0/VZ 294 : C 295 : C- Calcul de la vitesse de Chapman-Jouguet 296 : C 297 : VCJ = V1*SQRT((PAD-P1)/(V1-VAD)) 298 : C 299 : IF (RAD.GT.0.D0) IRET = 0 300 : RETURN 301 : END 302 : 303 : 304 :
© Cast3M 2003 - All rights reserved.
Disclaimer