Esope source files

   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 :