Télécharger la plaquette Cast3M


fr en

Télécharger xfpa.eso

Retour à la liste

Numérotation des lignes : 

xfpa

  1. C XFPA      SOURCE    CHAT      05/01/13    04:13:32     5004
  2.       SUBROUTINE XFPA(ANU,YP,ROG,RAP,N,UET,VNORM,AK)
  3. C*****************************************************************
  4. C
  5. C OBJET : calcul de la vitesse de depot d'aerosols AK dans chaque
  6. C         element de la ligne paroi.
  7. C
  8. C SYNTAXE : AK = FPA ANU YP UET VNORM ROG RAP
  9. C
  10. C AK   : CHPOINT SCAL CENTRE (vitesse de depot)
  11. C ANU  : FLOTTANT (viscosite)
  12. C YP   : FLOTTANT (epaisseur de la couche limite)
  13. C UET  : CHPOINT SCAL CENTRE (vitesse de frottement)
  14. C VNORM : CHPOINT VECT FACE (normale a la paroi)
  15. C ROG  : POINT (masse volumique * g)
  16. C RAP  : FLOTTANT (rayon des particules)
  17. C
  18. C EST APPELEE PAR KFPA.SOR
  19. C*****************************************************************
  20. C
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8 (A-H,O-Z)
  23.       DIMENSION UET(*),VNORM(*),AK(*),ROG(*)
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. C
  28. C
  29. C
  30. C
  31. C ---------- parametres de la particule ------------------------------
  32. C - coefficient de correction de Cunningham :
  33.       CUN = 1.0D0 + 8.296D-8/RAP + 2.64D-8/(RAP*EXP(1.667D7*RAP))
  34. C
  35. C - masse volumique ROP :
  36.       ROP = ROG(1)*ROG(1)
  37.       DO 20 I=2,IDIM
  38.             ROP = ROP + ROG(I)*ROG(I)
  39.  20   CONTINUE
  40.       ROP = SQRT(ROP)/9.81D0
  41. C
  42. C - vitesse de sedimentation incomplete VSI :
  43.       VSI = 1.23D4*RAP*RAP*CUN
  44.       VS  = VSI*ROP*9.81D0
  45. C
  46. C - temps de relaxation TAU :
  47.       TAU = VSI*ROP
  48. C
  49. C - nombre de SCHMIDT :
  50.       SCHMIDT = 5.62D21*ANU*ANU*RAP/CUN
  51. C
  52. C - constante de la vitesse de SAFFMAN (adimensionnelle) :
  53.       CTSAF = 2.23D5*RAP*TAU
  54. C
  55.             DO 3453 IEL=1,N
  56. C
  57. C------- calcul de AKS = vitesse de sedimentation --------------------
  58. C          ( d' apres l'orientation de l'element )
  59. C
  60.       PROSCA = ROG(1)*VNORM(IEL)
  61.       DO 30 I=2,IDIM
  62.             PROSCA = PROSCA + ROG(I)*VNORM(IEL+(I-1)*N)
  63.  30   CONTINUE
  64.       IF (PROSCA.LT.0.0D0) PROSCA = 0.0D0
  65.       AKS = VSI*PROSCA
  66. C
  67. C
  68. C------- calcul de AKDI = vitesse de depot par diffusion + inertie ----
  69. C
  70.       UIT = UET(IEL)
  71.       RAPLUS = RAP*UIT/ANU
  72.       TOPLUS = TAU*UIT*UIT/ANU
  73.       YFIN = YP*UIT/ANU
  74. C
  75. C - initialisations et boucle :
  76.       YPLUS1 = RAPLUS
  77.       YPLUS2 = RAPLUS
  78.       YPAS = 1.0D-4
  79.       ANUT = 0.0D0
  80.       ZETA1 = SCHMIDT
  81.       ZETA2 = SCHMIDT
  82.       ZETA = SCHMIDT*1.0D-4
  83.       VPRIMF = 0.0D0
  84.       VPRIMP1 = 0.0D0
  85.       VPRIMP2 = 0.0D0
  86.       TLPLUS = 0.0D0
  87.       VMIGT = 0.0D0
  88.       VSAFF = CTSAF
  89.       VPLUS = 0.0D0
  90.       CPLUS = 0.0D0
  91. C
  92. C
  93.       DO 10   I = 1,200
  94.               YPLUS2 = YPLUS2 + YPAS
  95.          IF (YPLUS2.LE.YFIN) THEN
  96.              IF (YPLUS2.LT.3) THEN
  97.                  ANUT = (YPLUS2/1.115D1)**3
  98.                  VSAFF = CTSAF/SQRT(1.0D0 + (YPLUS2/1.14D1)**2)
  99.                ELSEIF (YPLUS2.LT.52) THEN
  100.                  ANUT = (YPLUS2/1.14D1)**2 - 4.9774D-2
  101.                  VSAFF = CTSAF/SQRT(1.0D0 + (YPLUS2/1.14D1)**2)
  102.                ELSE
  103.                  ANUT = 4.0D-1*YPLUS2
  104.                  VSAFF = CTSAF/SQRT(1.0D0 + 4.0D-1*YPLUS2)
  105.              ENDIF
  106. C
  107.              IF (YPLUS2.LT.30) THEN
  108.                  VPRIM = 3.3D-2*(1.0D0 - EXP(-YPLUS2/3.837D0))
  109.                  VPRIM = VPRIM*YPLUS2
  110.                  VPR30 = VPRIM
  111.                ELSE
  112.                  IF (YFIN.EQ.3.0D1)   YFIN = 3.1D1
  113.                  YPLUS3 = (YPLUS2-3.0D1)/(YFIN-3.0D1)
  114.                  VPRIM = VPR30 - (VPR30 - 1.82D0)*YPLUS3
  115.              ENDIF
  116.                  VPRIMF = VPRIM*VPRIM
  117. C
  118.              ZETA2 = SCHMIDT/(1.0D0 + ANUT*SCHMIDT)
  119.              ZETA = (ZETA1 + ZETA2)*YPAS/2.0D0
  120.                    ZETA1 = ZETA2
  121. C
  122.              TLPLUS = ANUT/VPRIMF
  123.              VPRIMP2 = VPRIMF/(1.0D0 + TOPLUS/TLPLUS)
  124.              VMIGT = TOPLUS*(VPRIMP2 - VPRIMP1)/YPAS
  125.                    VPRIMP1 = VPRIMP2
  126.              VPLUS = VMIGT + VSAFF
  127.              CPLUS = CPLUS + ZETA*(1.0D0 - CPLUS*VPLUS)
  128. C
  129. C       WRITE(6,*) 'CPLUS =' , CPLUS
  130. C       WRITE(6,*) 'ZETA =' , ZETA
  131. C       WRITE(6,*) 'VPLUS =' , VPLUS
  132. C       WRITE(6,*) 'VPRIMP =' , VPRIMP2
  133. C       WRITE(6,*) 'YPLUS =' , YPLUS2
  134.              YPAS = (YPLUS2 - YPLUS1)*1.1D0
  135.                    YPLUS1 = YPLUS2
  136.          ENDIF
  137.   10  CONTINUE
  138. C
  139.       AKPLUS = 1.0D0/CPLUS
  140.       AKDI = AKPLUS*UIT
  141. C
  142. C
  143. C------- calcul de AKT = vitesse de depot totale ----------------------
  144. C
  145.       AKT = AKS + AKDI
  146.       AK(IEL) = AKT
  147. C
  148. C
  149.  3453 CONTINUE
  150. C
  151.       RETURN
  152.       END
  153.  
  154.  

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