Télécharger la plaquette Cast3M


fr en

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