Télécharger deto2.eso

Retour à la liste

Numérotation des lignes :

  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
  105. CALL DETOH(T1,X,AN,AH2,AH2O,AN2,AO2,TINT,H1,IINT,NBC,NBT)
  106. CALL DETORO(P1,T1,X,AN,R1,NBC)
  107. V1 = 1.D0/R1
  108. C
  109. C- On verifie que le mélange est susceptible de déroner.
  110. C- Dans le cas contraire, on renvoie la condition initiale
  111. C
  112. IF ((X(1).LE.0.04D0) .OR. (X(3).GE.0.6D0)) THEN
  113. RAD = R1
  114. TAD = T1
  115. PAD = P1
  116. FACT = 0.D0
  117. VCJ = 0.D0
  118. RZ = R1
  119. TZ = T1
  120. PZ = P1
  121. RAICC = R1
  122. TAICC = T1
  123. PAICC = P1
  124. IRET = 0
  125. RETURN
  126. ENDIF
  127. C
  128. C- Initialisation du taux de combustion (1=complète)
  129. C
  130. FACT = 1.D0
  131. AICC = .TRUE.
  132. C---------------------------------------------------------------------
  133. C- Retour du NEWTON sur le taux de combustion : la variable pivot est
  134. C- le taux de combustion.
  135. C---------------------------------------------------------------------
  136. 70 CONTINUE
  137. C
  138. C- Evaluation des fractions molaires pour un taux de combustion égal à
  139. C- FACT et prise en compte d'un éventuel déficit en oxygène du mélange
  140. C
  141. XC(1) = (X(1) - FACT*X(1)) * P1/T1
  142. XC(2) = (X(2) - 0.5D0*FACT*X(1)) * P1/T1
  143. XC(3) = (X(3) + FACT*X(1)) * P1/T1
  144. XC(4) = X(4) * P1/T1
  145. IF (XC(2).LT.0.D0) THEN
  146. XC(2) = 0.D0
  147. XC(1) = (X(1) - 2.D0*X(2)) * P1/T1
  148. XC(3) = (X(3) + 2.D0*X(2)) * P1/T1
  149. ENDIF
  150. XTOT = 0.D0
  151. DO 3 K = 1,NBC
  152. XTOT = XTOT + XC(K)
  153. 3 CONTINUE
  154. DO 4 K = 1,NBC
  155. XC(K) = XC(K) / XTOT
  156. 4 CONTINUE
  157. C
  158. C- Initialisation de la densité R0 à la densité initiale
  159. C
  160. R0 = R1
  161. C---------------------------------------------------------------------
  162. C- Retour du NEWTON sur le test d'appartenance du point initial à
  163. C- la tangente à l'hugoniot passant par le point CJ candidat :
  164. C- la variable pivot est la densité R0.
  165. C---------------------------------------------------------------------
  166. 60 CONTINUE
  167. R2 = R0
  168. INIT=.TRUE.
  169. DROITE = .FALSE.
  170. GAUCHE = .FALSE.
  171. C---------------------------------------------------------------------
  172. C- Retour afin de calculer les états à gauche et à droite du point
  173. C- (PC,TC,R2) susceptible d'etre l'etat de CJ afin d'évaluer la
  174. C- tangente à l'hugoniot.
  175. C---------------------------------------------------------------------
  176. 50 CONTINUE
  177. PC = P1
  178. C---------------------------------------------------------------------
  179. C- Retour du NEWTON permettant d'évaluer l'état d'équilibre (PC,TC)
  180. C- la densité étant fixée (=R2=R0) : on cherche (PC,TC) appartenant
  181. C- à la courbe d'hugoniot pour un taux de réaction égal à FACT. La
  182. C- variable pivot est la pression PC (de facon indireste TC).
  183. C---------------------------------------------------------------------
  184. 20 CONTINUE
  185. WTM2 = XC(1)*AN(1) + XC(2)*AN(2) + XC(3)*AN(3) + XC(4)*AN(4)
  186. WTM2 = WTM2 * 1.D-3
  187. TC = WTM2*PC/(R2*8.314D0)
  188. CALL DETOH(TC,XC,AN,AH2,AH2O,AN2,AO2,TINT,H2,IINT,NBC,NBT)
  189. FX = 2.D0*(H2-H1)-((PC-P1)*((1.D0/R1)+(1.D0/R2)))
  190. X1 = 1.D0
  191. X2 = 5.D-5
  192. CALL NEWTON(PC,FX,X1,X2,1000,3,NFIN3,NIT,NITER,TCX1,TCY1,TCX2,
  193. $ TCY2)
  194. IF (NFIN3 .EQ. 0) GOTO 20
  195. IF (NFIN3 .EQ. 2) RETURN
  196. C
  197. C- Le point (PC,TC,R2) trouvé est l'état CJ à condition que la pente
  198. C- à l'hugoniot passe par le point initial (P1,T1,R1). Pour évaluer la
  199. C- pente, on calcule un état gauche et un état droit de (PC,TC,R2) en
  200. C- modifiant la densité de 1% (par défaut et par excès)
  201. C
  202. IF (GAUCHE) THEN
  203. PG = PC
  204. VG = 1.D0/R2
  205. GAUCHE = .FALSE.
  206. ENDIF
  207. IF (DROITE) THEN
  208. PD = PC
  209. VD = 1.D0/R2
  210. DROITE = .FALSE.
  211. GAUCHE = .TRUE.
  212. ENDIF
  213. IF (INIT) THEN
  214. PAD = PC
  215. TAD = TC
  216. VAD = 1.D0/R2
  217. C
  218. C- La première convergence donne l'état AICC puisque le point obtenu
  219. C- est le point sur l'hugoniot correspondant à un taux de combustion
  220. C- complète et à la densité initiale.
  221. C
  222. IF (AICC) THEN
  223. PAICC = PC
  224. TAICC = TC
  225. VAICC = 1.D0/R2
  226. AICC = .FALSE.
  227. ENDIF
  228. INIT = .FALSE.
  229. DROITE = .TRUE.
  230. ENDIF
  231. IF (DROITE) THEN
  232. R2 = R0*0.99D0
  233. GOTO 50
  234. ENDIF
  235. IF (GAUCHE) THEN
  236. R2 = R0*1.01D0
  237. GOTO 50
  238. ENDIF
  239. C
  240. C- Equation de la droite passant par les points (PTEST,V1) et (PC,VC)
  241. C- P = ALP*V + BCS
  242. C- On vérifie l'ecart entre PTEST et P1 ; en cas de non convergence
  243. C- on définit un nouveau R0 et on recommence la détermination d'un
  244. C- état CJ potentiel.
  245. C
  246. ALP = (PD-PG)/(VD-VG)
  247. BCS = PAD - (ALP*VAD)
  248. PTEST = ALP*V1+BCS
  249. FY = PTEST - P1
  250. X11 = 1.D0
  251. X21 = 5.D-5
  252. CALL NEWTON(R0,FY,X11,X21,1000,2,NFIN2,NIT,NITER,TCX1,TCY1,
  253. $ tcx2,tcy2)
  254. IF (NFIN2 .EQ. 0) GOTO 60
  255. IF (NFIN2 .EQ. 2) RETURN
  256. C
  257. C- XC, PAD, TAD et VAD caractérisent un état CJ qui est le bon à
  258. C- condition que l'équilibre chimique soit réalisé : On calcule
  259. C- le taux de production de H2 (i.e. de chaque espece puisque ici
  260. C- cinétique à une réaction) qui doit etre nul ; sinon, on diminue
  261. C- le taux de combustion.
  262. C
  263. CALL DETOWP(PAD,TAD,XC,AH2,AH2O,AO2,TINT,WDOT,IINT,NBC,NBT)
  264. FZ = WDOT(1)
  265. X13 =-0.001D0
  266. X23 = 1.D-20
  267. CALL NEWTON(FACT,FZ,X13,X23,1000,1,NFIN1,NIT,NITER,TCX1,TCY1,
  268. $ TCX2,TCY2)
  269. IF (NFIN1 .EQ. 0) GOTO 70
  270. IF (NFIN1 .EQ. 2) RETURN
  271. C
  272. C- Initialisation du volume spécifique de l'état ZND et de la pente
  273. C- de la droite de Rayleigh
  274. C
  275. VZ = V1/5.D0
  276. PENTE = (PAD-P1) / (VAD-V1)
  277. C---------------------------------------------------------------------
  278. C- Retour du NEWTON du calcul ZND : on cherche le point d'intersection
  279. C- entre l'hugoniot passant par le point initial et la droite passant
  280. C- part CJ et le point initial. La variable pivot est le volume
  281. C- spécifique VZ.
  282. C---------------------------------------------------------------------
  283. 80 CONTINUE
  284. PZ = P1 + (PENTE*(VZ-V1))
  285. WTM1 = X(1)*AN(1) + X(2)*AN(2) + X(3)*AN(3) + X(4)*AN(4)
  286. WTM1 = WTM1 * 1.D-3
  287. TZ = WTM1*PZ*VZ / 8.314D0
  288. CALL DETOH(TZ,X,AN,AH2,AH2O,AN2,AO2,TINT,HZ,IINT,NBC,NBT)
  289. FQ =-2.D0 * (HZ-H1) + (PZ-P1)*(V1+VZ)
  290. X1 =-0.01D0
  291. X2 = 5.D-5
  292. X2 = 5.D-10
  293. CALL NEWTON(VZ,FQ,X1,X2,1000,1,NFIN,NIT,NITER,TCX1,TCY1,
  294. $ TCX2,TCY2)
  295. IF (NFIN .EQ. 0) GOTO 80
  296. IF (NFIN .EQ. 2) RETURN
  297. C
  298. C- Calcul des densités à partir des volumes spécifiques
  299. C
  300. RAD = 1.D0/VAD
  301. RAICC = 1.D0/VAICC
  302. RZ = 1.D0/VZ
  303. C
  304. C- Calcul de la vitesse de Chapman-Jouguet
  305. C
  306. VCJ = V1*SQRT((PAD-P1)/(V1-VAD))
  307. C
  308. IF (RAD.GT.0.D0) IRET = 0
  309. RETURN
  310. END
  311.  
  312.  
  313.  
  314.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales