Télécharger flacr3.eso

Retour à la liste

Numérotation des lignes :

flacr3
  1. C FLACR3 SOURCE OF166741 24/12/13 21:15:52 12097
  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 CSI<EPSCSI ou CSI>1+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.  
  75. -INC PPARAM
  76. -INC CCOPTIO
  77. -INC SMLREEL
  78. POINTEUR MLRECO.MLREEL, MLRMAS.MLREEL, MLRH0K.MLREEL
  79. -INC SMELEME
  80. POINTEUR MELCEN.MELEME, MELEFL.MELEME
  81. -INC SMCHPOI
  82. INTEGER N, NC
  83. POINTEUR MPOVRO.MPOVAL, MPOVAY.MPOVAL,MPOVYI.MPOVAL,
  84. & MPOVYF.MPOVAL, MPOVC.MPOVAL, MPODRE.MPOVAL, MPODRY.MPOVAL,
  85. & MPOCSI.MPOVAL,MPOCRI.MPOVAL,MPOVDX.MPOVAL
  86. -INC SMLENTI
  87. POINTEUR MLECEN.MLENTI
  88. C
  89. C**** Les variables
  90. C
  91. INTEGER MCEN,MFACEL,IRC,IYC,IYINIT,IYFINA,IVCAR,ICHRET,ICHRY
  92. $ ,IGEOM,ICEN,NCEN,NFAC,NLCF,NGCEG,NGCED,NLCEG,NLCED
  93. & ,IESP,NESP,IDX
  94. REAL*8 EPSI, DELTAT, Y1F, Y1I, Y1, VCSIG, VCSID, VCSI2G, VCSI2D
  95. & ,EPS12, DY, DYMAX, DCSI, DCSI1, DY1, RHO, DY2, YMAX, DYMAX1
  96. & ,ERRTOL, EPSCSI, CSIG
  97. CHARACTER*8 TYPE
  98. PARAMETER(ERRTOL=1.0D-6)
  99. C
  100. NESP=MLRH0K.PROG(/1)
  101. MELCEN=MCEN
  102. MELEFL=MFACEL
  103. C
  104. C**** KRIPAD pour la correspondance global/local de centre
  105. C
  106. CALL KRIPAD(MELCEN,MLECEN)
  107. IF(IERR .NE. 0)GOTO 9999
  108. C SEGINI MLECEN
  109. SEGACT MELCEN
  110. NCEN=MELCEN.NUM(/2)
  111. SEGACT MELEFL
  112. NFAC=MELEFL.NUM(/2)
  113. C
  114. C**** Lectures de CHPOINT
  115. C
  116. CALL LICHT(IRC,MPOVRO,TYPE,IGEOM)
  117. C SEGACT MPOVRO*MOD
  118. CALL LICHT(IYC,MPOVAY,TYPE,IGEOM)
  119. C SEGACT MPOVAY*MOD
  120. CALL LICHT(IYINIT,MPOVYI,TYPE,IGEOM)
  121. C SEGACT MPOVYI*MOD
  122. CALL LICHT(IYFINA,MPOVYF,TYPE,IGEOM)
  123. C SEGACT MPOVYF*MOD
  124. CALL LICHT(IVCAR,MPOVC,TYPE,IGEOM)
  125. C SEGACT MPOVYF*MOD
  126. CALL LICHT(IDX,MPOVDX,TYPE,IGEOM)
  127. C SEGACT MPOVDX*MOD
  128. CALL LICHT(ICHRET,MPODRE,TYPE,IGEOM)
  129. C SEGACT MPODRE*MOD
  130. CALL LICHT(ICHRY,MPODRY,TYPE,IGEOM)
  131. C SEGACT MPODRY*MOD
  132. C
  133. C**** Creation du MPOVAL qui contient csi
  134. C Creation du MPOVAL du critere
  135. C
  136. N=NCEN
  137. NC=1
  138. SEGINI MPOCSI
  139. SEGINI MPOCRI
  140. C
  141. C**** Calcul de DYMAX pour la premiere espece
  142. C Controle de la densité positive
  143. C Controle des fractions massiques (Yi>0, sum_i Y_i < 1)
  144. C
  145. DYMAX = 0.0D0
  146. YMAX=0.0D0
  147. DO ICEN=1,NCEN,1
  148. Y1F = MPOVYF.VPOCHA(ICEN,1)
  149. Y1I = MPOVYI.VPOCHA(ICEN,1)
  150. DY = Y1F - Y1I
  151. DYMAX=MAX(ABS(DY),DYMAX)
  152. Y1=MPOVAY.VPOCHA(ICEN,1)
  153. IF(Y1 .LT. 0.0D0)THEN
  154. WRITE(IOIMP,*) 'CHPO2 < 0 ???'
  155. CALL ERREUR(21)
  156. GOTO 9999
  157. ELSE
  158. YMAX=Y1
  159. ENDIF
  160. RHO=MPOVRO.VPOCHA(ICEN,1)
  161. IF(RHO .LT. 0)THEN
  162. WRITE(IOIMP,*) 'CHPO1 < 0 ???'
  163. CALL ERREUR(21)
  164. GOTO 9999
  165. ENDIF
  166. RHO=MPOVYI.VPOCHA(ICEN,1)
  167. IF((RHO .LT. 0) .OR. (RHO .GT. 1))THEN
  168. 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???'
  175. CALL ERREUR(21)
  176. GOTO 9999
  177. ENDIF
  178. RHO=MPOVC.VPOCHA(ICEN,1)
  179. IF((RHO .LT. 0))THEN
  180. WRITE(IOIMP,*) 'CHPO5 < 0'
  181. CALL ERREUR(21)
  182. GOTO 9999
  183. ENDIF
  184. DO IESP=2,NESP,1
  185. Y1=MPOVAY.VPOCHA(ICEN,IESP)
  186. IF(Y1 .LT. 0.0D0)THEN
  187. WRITE(IOIMP,*) 'CHPO2 < 0 ???'
  188. CALL ERREUR(21)
  189. GOTO 9999
  190. ELSE
  191. YMAX=YMAX+Y1
  192. ENDIF
  193. ENDDO
  194. IF(YMAX .GT. 1.0D0)THEN
  195. WRITE(IOIMP,*) 'sum CHPO2 > 1 ???'
  196. CALL ERREUR(21)
  197. GOTO 9999
  198. ENDIF
  199. ENDDO
  200. C
  201. C**** Calcul de CSI \in (0,1)
  202. C
  203. DO ICEN=1,NCEN,1
  204. Y1F = MPOVYF.VPOCHA(ICEN,1)
  205. Y1I = MPOVYI.VPOCHA(ICEN,1)
  206. DY = Y1F - Y1I
  207. IF(ABS(DY) .LE. (ERRTOL*DYMAX))THEN
  208. C
  209. C********** Pas de combustion en cette region
  210. C
  211. MPOCSI.VPOCHA(ICEN,1)=0.0D0
  212. C
  213. ELSE
  214. Y1 = MPOVAY.VPOCHA(ICEN,1)
  215. VCSIG = (Y1 - Y1I) / DY
  216. C
  217. C********** On n'accepte pas csi > 1.1 or csi < -0.1
  218. C
  219. IF((VCSIG .GT. (1.0D0+EPSCSI)) .OR.
  220. & (VCSIG .LT. (-1*EPSCSI)))THEN
  221. WRITE(IOIMP,*) 'Progress variable = ???'
  222. C 21 2
  223. C Données incompatibles
  224. CALL ERREUR(21)
  225. GOTO 9999
  226. ELSEIF(VCSIG .GT. 1.0D0)THEN
  227. VCSIG = 1.0D0
  228. ELSEIF(VCSIG .LT. 0.0D0)THEN
  229. VCSIG = 0.0D0
  230. ENDIF
  231. MPOCSI.VPOCHA(ICEN,1)= VCSIG
  232. ENDIF
  233. ENDDO
  234. C
  235. C**** Le critere CREBCOM
  236. C
  237. DO NLCF = 1, NFAC, 1
  238. C
  239. C******* NLCF = numero local du centre de facel
  240. C NGCEG = numero global du centre ELT "gauche"
  241. C NLCEG = numero local du centre ELT "gauche"
  242. C NGCED = numero global du centre ELT "droite"
  243. C NLCED = numero local du centre ELT "droite"
  244. C
  245. NGCEG = MELEFL.NUM(1,NLCF)
  246. NGCED = MELEFL.NUM(3,NLCF)
  247. NLCEG = MLECEN.LECT(NGCEG)
  248. NLCED = MLECEN.LECT(NGCED)
  249. C
  250. VCSIG=MPOCSI.VPOCHA(NLCEG,1)
  251. VCSID=MPOCSI.VPOCHA(NLCED,1)
  252. VCSI2G=VCSIG*VCSIG
  253. VCSI2D=VCSID*VCSID
  254. C
  255. IF(NLCEG .EQ. NLCED)THEN
  256. C
  257. C********** Murs
  258. C
  259. MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) + (0.5D0 *
  260. & VCSI2D)
  261. C
  262. ELSE
  263. C
  264. MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) +
  265. & (VCSI2D - (0.5D0 * VCSI2G))
  266. MPOCRI.VPOCHA(NLCED,1)=MPOCRI.VPOCHA(NLCED,1) +
  267. & (VCSI2G - (0.5D0 * VCSI2D))
  268. C
  269. ENDIF
  270. ENDDO
  271. C
  272. C**** Calcul des increments de l'energie et des densités
  273. C
  274. EPS12 = EPSI * EPSI
  275. DO ICEN = 1, NCEN, 1
  276. CSIG = MPOCSI.VPOCHA(ICEN,1)
  277. VCSIG = MPOCRI.VPOCHA(ICEN,1)
  278. C
  279. C******* In 2D, contribution of the ideal upper and lower cells
  280. C
  281. IF(IDIM .EQ. 2) VCSIG = VCSIG + (CSIG * CSIG)
  282. IF(VCSIG .GE. (EPS12*(1.0D0+ERRTOL)))THEN
  283. C
  284. C********** Il y a combustion
  285. C
  286. DCSI=MPOVC.VPOCHA(ICEN,1)*DELTAT/MPOVDX.VPOCHA(ICEN,1)
  287. DCSI1=1.0D0-CSIG
  288. DCSI=MIN(DCSI,DCSI1)
  289. DY1=(MPOVYF.VPOCHA(ICEN,1)-MPOVYI.VPOCHA(ICEN,1))*DCSI
  290. C
  291. C********** On force la positivité de MPOVAY.VPOCHA(ICEN,1)
  292. C
  293. IF(DY1 .GT. 0.0D0)THEN
  294. DY2 = 1.0D0 - MPOVAY.VPOCHA(ICEN,1)
  295. ELSE
  296. DY2 = MPOVAY.VPOCHA(ICEN,1)
  297. ENDIF
  298. DY1=SIGN(MIN(ABS(DY2),ABS(DY1)),DY1)
  299. C
  300. RHO=MPOVRO.VPOCHA(ICEN,1)
  301. MPODRY.VPOCHA(ICEN,1)=DY1*RHO
  302. MPODRE.VPOCHA(ICEN,1)=-1.0D0*DY1*RHO*MLRH0K.PROG(1)
  303. DY1=DY1/(MLRMAS.PROG(1)*MLRECO.PROG(1))
  304. DO IESP=2,NESP,1
  305. DY=DY1*(MLRMAS.PROG(IESP)*MLRECO.PROG(IESP))
  306. DYMAX1=ABS(DYMAX/(MLRMAS.PROG(1)*MLRECO.PROG(1))*
  307. & (MLRMAS.PROG(IESP)*MLRECO.PROG(IESP)))
  308. IF(DY .GT. 0)THEN
  309. DY2=1.0D0 - MPOVAY.VPOCHA(ICEN,IESP)
  310. ELSE
  311. DY2=MPOVAY.VPOCHA(ICEN,IESP)
  312. ENDIF
  313. C
  314. C************* On force la positivité de MPOVAY.VPOCHA(ICEN,IESP)
  315. C
  316. IF((ABS(DY)-ABS(DY2)).GE.(0.5D0*DYMAX1))THEN
  317. WRITE(IOIMP,*) 'CHPO2, CHPO3, CHPO4 = ???'
  318. C Données incompatibles
  319. CALL ERREUR(21)
  320. GOTO 9999
  321. ELSE
  322. DY=SIGN(MIN(ABS(DY2),ABS(DY)),DY)
  323. ENDIF
  324. MPODRY.VPOCHA(ICEN,IESP)=DY*RHO
  325. MPODRE.VPOCHA(ICEN,1)=MPODRE.VPOCHA(ICEN,1)-(DY*RHO
  326. $ *MLRH0K.PROG(IESP))
  327. ENDDO
  328. ENDIF
  329. ENDDO
  330. C
  331. SEGDES MELCEN
  332. SEGDES MELEFL
  333. SEGSUP MLECEN
  334. C
  335. SEGDES MPOVRO
  336. SEGDES MPOVAY
  337. SEGDES MPOVYI
  338. SEGDES MPOVYF
  339. SEGDES MPOVC
  340. SEGDES MPOVDX
  341. SEGDES MPODRE
  342. SEGDES MPODRY
  343. SEGSUP MPOCSI
  344. SEGSUP MPOCRI
  345. C
  346. 9999 RETURN
  347. END
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  

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