Télécharger xfpa.eso

Retour à la liste

Numérotation des lignes :

  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. -INC CCOPTIO
  25. C
  26. C
  27. C
  28. C
  29. C ---------- parametres de la particule ------------------------------
  30. C - coefficient de correction de Cunningham :
  31. CUN = 1.0D0 + 8.296D-8/RAP + 2.64D-8/(RAP*EXP(1.667D7*RAP))
  32. C
  33. C - masse volumique ROP :
  34. ROP = ROG(1)*ROG(1)
  35. DO 20 I=2,IDIM
  36. ROP = ROP + ROG(I)*ROG(I)
  37. 20 CONTINUE
  38. ROP = SQRT(ROP)/9.81D0
  39. C
  40. C - vitesse de sedimentation incomplete VSI :
  41. VSI = 1.23D4*RAP*RAP*CUN
  42. VS = VSI*ROP*9.81D0
  43. C
  44. C - temps de relaxation TAU :
  45. TAU = VSI*ROP
  46. C
  47. C - nombre de SCHMIDT :
  48. SCHMIDT = 5.62D21*ANU*ANU*RAP/CUN
  49. C
  50. C - constante de la vitesse de SAFFMAN (adimensionnelle) :
  51. CTSAF = 2.23D5*RAP*TAU
  52. C
  53. DO 3453 IEL=1,N
  54. C
  55. C------- calcul de AKS = vitesse de sedimentation --------------------
  56. C ( d' apres l'orientation de l'element )
  57. C
  58. PROSCA = ROG(1)*VNORM(IEL)
  59. DO 30 I=2,IDIM
  60. PROSCA = PROSCA + ROG(I)*VNORM(IEL+(I-1)*N)
  61. 30 CONTINUE
  62. IF (PROSCA.LT.0.0D0) PROSCA = 0.0D0
  63. AKS = VSI*PROSCA
  64. C
  65. C
  66. C------- calcul de AKDI = vitesse de depot par diffusion + inertie ----
  67. C
  68. UIT = UET(IEL)
  69. RAPLUS = RAP*UIT/ANU
  70. TOPLUS = TAU*UIT*UIT/ANU
  71. YFIN = YP*UIT/ANU
  72. C
  73. C - initialisations et boucle :
  74. YPLUS1 = RAPLUS
  75. YPLUS2 = RAPLUS
  76. YPAS = 1.0D-4
  77. ANUT = 0.0D0
  78. ZETA1 = SCHMIDT
  79. ZETA2 = SCHMIDT
  80. ZETA = SCHMIDT*1.0D-4
  81. VPRIMF = 0.0D0
  82. VPRIMP1 = 0.0D0
  83. VPRIMP2 = 0.0D0
  84. TLPLUS = 0.0D0
  85. VMIGT = 0.0D0
  86. VSAFF = CTSAF
  87. VPLUS = 0.0D0
  88. CPLUS = 0.0D0
  89. C
  90. C
  91. DO 10 I = 1,200
  92. YPLUS2 = YPLUS2 + YPAS
  93. IF (YPLUS2.LE.YFIN) THEN
  94. IF (YPLUS2.LT.3) THEN
  95. ANUT = (YPLUS2/1.115D1)**3
  96. VSAFF = CTSAF/SQRT(1.0D0 + (YPLUS2/1.14D1)**2)
  97. ELSEIF (YPLUS2.LT.52) THEN
  98. ANUT = (YPLUS2/1.14D1)**2 - 4.9774D-2
  99. VSAFF = CTSAF/SQRT(1.0D0 + (YPLUS2/1.14D1)**2)
  100. ELSE
  101. ANUT = 4.0D-1*YPLUS2
  102. VSAFF = CTSAF/SQRT(1.0D0 + 4.0D-1*YPLUS2)
  103. ENDIF
  104. C
  105. IF (YPLUS2.LT.30) THEN
  106. VPRIM = 3.3D-2*(1.0D0 - EXP(-YPLUS2/3.837D0))
  107. VPRIM = VPRIM*YPLUS2
  108. VPR30 = VPRIM
  109. ELSE
  110. IF (YFIN.EQ.3.0D1) YFIN = 3.1D1
  111. YPLUS3 = (YPLUS2-3.0D1)/(YFIN-3.0D1)
  112. VPRIM = VPR30 - (VPR30 - 1.82D0)*YPLUS3
  113. ENDIF
  114. VPRIMF = VPRIM*VPRIM
  115. C
  116. ZETA2 = SCHMIDT/(1.0D0 + ANUT*SCHMIDT)
  117. ZETA = (ZETA1 + ZETA2)*YPAS/2.0D0
  118. ZETA1 = ZETA2
  119. C
  120. TLPLUS = ANUT/VPRIMF
  121. VPRIMP2 = VPRIMF/(1.0D0 + TOPLUS/TLPLUS)
  122. VMIGT = TOPLUS*(VPRIMP2 - VPRIMP1)/YPAS
  123. VPRIMP1 = VPRIMP2
  124. VPLUS = VMIGT + VSAFF
  125. CPLUS = CPLUS + ZETA*(1.0D0 - CPLUS*VPLUS)
  126. C
  127. C WRITE(6,*) 'CPLUS =' , CPLUS
  128. C WRITE(6,*) 'ZETA =' , ZETA
  129. C WRITE(6,*) 'VPLUS =' , VPLUS
  130. C WRITE(6,*) 'VPRIMP =' , VPRIMP2
  131. C WRITE(6,*) 'YPLUS =' , YPLUS2
  132. YPAS = (YPLUS2 - YPLUS1)*1.1D0
  133. YPLUS1 = YPLUS2
  134. ENDIF
  135. 10 CONTINUE
  136. C
  137. AKPLUS = 1.0D0/CPLUS
  138. AKDI = AKPLUS*UIT
  139. C
  140. C
  141. C------- calcul de AKT = vitesse de depot totale ----------------------
  142. C
  143. AKT = AKS + AKDI
  144. AK(IEL) = AKT
  145. C
  146. C
  147. 3453 CONTINUE
  148. C
  149. RETURN
  150. END
  151.  
  152.  

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