part.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

   ! Loop over all particle/droplet classes that have the given evaporative index

   PART_CLASS_LOOP: DO N_PC = 1,N_PART

      PC => PARTICLE_CLASS(N_PC)
      IF (PC%EVAP_INDEX/=EVAP_INDEX) CYCLE PART_CLASS_LOOP
      DROP_DEN = 0._EB
      DROP_TMP = 0._EB
      DROP_RAD = 0._EB
      C_DROP   = PC%C_P
      FTPR     = PC%FTPR
      TMP_MELT = PC%TMP_MELT
      TMP_BOIL = PC%TMP_V
      IGAS     = PC%SPECIES_INDEX
      MW_DROP  = SPECIES(IGAS)%MW
      MW_RATIO = SPECIES(0)%MW/MW_DROP
      H_V      = PC%H_V
      DHOR     = H_V*MW_DROP/R0

      ! Loop through all droplets in the class and determine the depth of the liquid film on each surface cell

      MASS_SUMMING_LOOP: DO I=1,NLP
         DR => DROPLET(I)
         IF (DR%IOR==0)        CYCLE MASS_SUMMING_LOOP
         IF (DR%WALL_INDEX==0) CYCLE MASS_SUMMING_LOOP
         IF (DR%CLASS /= N_PC) CYCLE MASS_SUMMING_LOOP
         IF (DR%R<=0._EB)      CYCLE MASS_SUMMING_LOOP
         IW = DR%WALL_INDEX
         FILM_THICKNESS(IW) = FILM_THICKNESS(IW) + DR%PWT*DR%R**3/AW(IW)
      ENDDO MASS_SUMMING_LOOP
      FILM_THICKNESS = FILM_THICKNESS*FTPR/PC%DENSITY

      FILM_THICKNESS = MAX(MINIMUM_FILM_THICKNESS,FILM_THICKNESS) 

      ! Loop through all droplets within the class and determine mass/energy transfer

      DROPLET_LOOP: DO I=1,NLP

         DR => DROPLET(I)
         IF (DR%CLASS /= N_PC) CYCLE DROPLET_LOOP
         IF (DR%R<=0._EB)      CYCLE DROPLET_LOOP

         ! Determine the current coordinates of the particle

         XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
         YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
         ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
         II  = FLOOR(XI+1._EB)
         JJ  = FLOOR(YJ+1._EB)
         KK  = FLOOR(ZK+1._EB)
         
         ! Initialize droplet thermophysical data

         R_DROP   = DR%R
         M_DROP   = FTPR*R_DROP**3
         TMP_DROP = DR%TMP
         WGT      = DR%PWT
               
         ! Gas conditions

         TMP_G  = TMP(II,JJ,KK)
         RHO_G  = RHO(II,JJ,KK)
         MU_AIR = SPECIES(0)%MU(MIN(500,NINT(0.1_EB*TMP_G)))
         RVC    = RDX(II)*RDY(JJ)*RDZ(KK)
         M_GAS  = RHO_G/RVC        
         K_AIR  = CPOPR*MU_AIR
         IF (IGAS>0) THEN
            Y_GAS = YY(II,JJ,KK,IGAS)
         ELSE
            Y_GAS = 0._EB
         ENDIF

         ! Set variables for heat transfer on solid

         SOLID_OR_GAS_PHASE: IF (DR%IOR/=0 .AND. DR%WALL_INDEX>0) THEN

            IW   = DR%WALL_INDEX
            A_DROP = M_DROP/(FILM_THICKNESS(IW)*PC%DENSITY)
            TMP_WALL = TMP_F(IW) 
            SELECT CASE(ABS(DR%IOR))
               CASE(1)
                  V2 = 0.25_EB*(V(II,JJ,KK)+V(II,JJ-1,KK))**2
                  W2 = 0.25_EB*(W(II,JJ,KK)+W(II,JJ,KK-1))**2
                  VEL = SQRT(V2+W2)
               CASE(2)
                  U2 = 0.25_EB*(U(II,JJ,KK)+U(II-1,JJ,KK))**2
                  W2 = 0.25_EB*(W(II,JJ,KK)+W(II,JJ,KK-1))**2
                  VEL = SQRT(U2+W2)
               CASE(3)
                  U2 = 0.25_EB*(U(II,JJ,KK)+U(II-1,JJ,KK))**2
                  V2 = 0.25_EB*(V(II,JJ,KK)+V(II,JJ-1,KK))**2
                  VEL = SQRT(U2+V2)
            END SELECT
            LENGTH   = 1._EB
            RE_L     = MAX(5.E5_EB,RHO_G*VEL*LENGTH/MU_AIR)
            NUSSELT  = NU_FAC_WALL*RE_L**0.8_EB
            SHERWOOD = SH_FAC_WALL*RE_L**0.8_EB
            H_HEAT   = NUSSELT*K_AIR/LENGTH
            IF (PC%EVAPORATE) THEN
               H_MASS = SHERWOOD*D_AIR/LENGTH
            ELSE
               H_MASS = 0._EB
            ENDIF
            H_WALL   = H_SOLID
!            Q_DOT_RAD = A_DROP*(QRADIN(IW)-QRADOUT(IW))
            Q_DOT_RAD = A_DROP*QRADIN(IW)

         ELSE SOLID_OR_GAS_PHASE

            A_DROP   = 4._EB*PI*R_DROP**2
            NUSSELT  = 2._EB + NU_FAC_GAS*SQRT(DR%RE)
            SHERWOOD = 2._EB + SH_FAC_GAS*SQRT(DR%RE)
            H_HEAT   = NUSSELT *K_AIR/(2._EB*R_DROP)
            IF (PC%EVAPORATE) THEN
               H_MASS = SHERWOOD*D_AIR/(2._EB*R_DROP)
            ELSE
               H_MASS = 0._EB
            ENDIF
            H_WALL   = 0._EB
            TMP_WALL = TMPA
            IF (AVG_DROP_DEN(II,JJ,KK,EVAP_INDEX )>0._EB) THEN
               Q_DOT_RAD = QR_W(II,JJ,KK)/SUM(AVG_DROP_DEN(II,JJ,KK,:))*M_DROP
            ELSE
               Q_DOT_RAD = 0._EB
            ENDIF

         ENDIF SOLID_OR_GAS_PHASE
         
         ! Compute equilibrium droplet vapor mass fraction, Y_DROP, and its derivative w.r.t. droplet temperature
    
         IF (PC%EVAPORATE) THEN
            X_DROP  = MIN(1._EB,EXP(DHOR*(1._EB/TMP_BOIL-1._EB/TMP_DROP)))
            Y_DROP  = X_DROP/(MW_RATIO + (1._EB-MW_RATIO)*X_DROP)
            IF (TMP_DROP < TMP_BOIL) THEN
               DY_DTMP_DROP = (MW_RATIO/(X_DROP*(1._EB-MW_RATIO)+MW_RATIO)**2)*DHOR*X_DROP/TMP_DROP**2
            ELSE
               DY_DTMP_DROP = 0._EB
            ENDIF
         ELSE
            X_DROP = 0._EB
            Y_DROP = 0._EB
            DY_DTMP_DROP = 0._EB
         ENDIF

         ! Update the droplet temperature semi_implicitly
    
         DENOM = 1._EB + (H_HEAT + H_WALL + H_MASS*RHO_G*H_V*DY_DTMP_DROP)*DT*A_DROP/(2._EB*M_DROP*C_DROP) 

         TMP_DROP_NEW = ( TMP_DROP + DT*( Q_DOT_RAD + &
                          A_DROP*(H_HEAT*(TMP_G   -0.5_EB*TMP_DROP) + H_WALL*(TMP_WALL-0.5_EB*TMP_DROP) -  &
                          H_MASS*RHO_G*H_V*(Y_DROP-0.5_EB*DY_DTMP_DROP*TMP_DROP-Y_GAS))/(M_DROP*C_DROP)) ) / DENOM

         ! Compute the total amount of heat extracted from the gas, wall and radiative fields

         Q_RAD      = DT*Q_DOT_RAD
         Q_CON_GAS  = DT*A_DROP*H_HEAT*(TMP_G   -0.5_EB*(TMP_DROP+TMP_DROP_NEW))
         Q_CON_WALL = DT*A_DROP*H_WALL*(TMP_WALL-0.5_EB*(TMP_DROP+TMP_DROP_NEW))

         ! Compute the total amount of liquid evaporated
  
         M_VAP  = DT*A_DROP*H_MASS*RHO_G*(Y_DROP+0.5_EB*DY_DTMP_DROP*(TMP_DROP_NEW-TMP_DROP)-Y_GAS) 
         M_VAP = MAX(0._EB,MIN(M_VAP,M_DROP))
         IF (PC%EVAPORATE .AND. R_DROP<0.5_EB*PC%MINIMUM_DIAMETER) THEN  ! Evaporate small droplets
            M_VAP      = M_DROP
            Q_TEMP = Q_RAD+Q_CON_GAS+Q_CON_WALL
            Q_TEMP = M_VAP*H_V/Q_TEMP
            Q_CON_GAS  = Q_CON_GAS*Q_TEMP
            Q_CON_WALL = Q_CON_WALL*Q_TEMP
            Q_RAD      = Q_RAD*Q_TEMP
         ENDIF

         ! If the droplet temperature drops below its freezing point, just reset it

         IF (PC%EVAPORATE .AND. TMP_DROP_NEW<TMP_MELT) TMP_DROP_NEW = TMP_MELT

         ! If the droplet temperature reaches boiling, use only enough energy from gas to vaporize liquid

         IF (PC%EVAPORATE .AND. TMP_DROP_NEW>=TMP_BOIL) THEN  
            Q_TEMP = Q_RAD+Q_CON_GAS+Q_CON_WALL
            M_VAP = MIN(M_DROP,M_VAP + (TMP_DROP_NEW - TMP_BOIL)*C_DROP*M_DROP/H_V)
            Q_TEMP = M_VAP*H_V/Q_TEMP
            TMP_DROP_NEW = TMP_BOIL
            Q_CON_GAS  = Q_CON_GAS*Q_TEMP
            Q_CON_WALL = Q_CON_WALL*Q_TEMP
            Q_RAD      = Q_RAD*Q_TEMP
         ENDIF

         ! Update droplet quantities
         
         M_DROP = M_DROP - M_VAP
         DR%R   = (M_DROP/FTPR)**ONTH
         DR%TMP = TMP_DROP_NEW

         ! Compute surface cooling and density

         IF (DR%IOR/=0 .AND. DR%WALL_INDEX>0) THEN
            WCPUA(IW,EVAP_INDEX) = WCPUA(IW,EVAP_INDEX) + OMRAF*WGT*RDT*(Q_RAD+Q_CON_WALL)*RAW(IW)
            WMPUA(IW,EVAP_INDEX) = WMPUA(IW,EVAP_INDEX) + OMRAF*WGT*M_DROP*RAW(IW)
         ENDIF
         
         ! Decrease temperature of the gas cell

         TMP(II,JJ,KK) = (M_GAS*CP_GAMMA*TMP_G - WGT*Q_CON_GAS)/(M_GAS*CP_GAMMA)  
         TMP(II,JJ,KK) = MIN(TMPMAX,MAX(TMPMIN,TMP(II,JJ,KK)))

         ! Compute contribution to the divergence

         D_VAP(II,JJ,KK) = D_VAP(II,JJ,KK) + WGT*RDT*RVC/RHO_G * ( M_VAP*MW_RATIO - Q_CON_GAS/(CP_GAMMA*TMP(II,JJ,KK)) )

         ! Adjust mass of evaporated liquid to account for different Heat of Combustion between fuel droplet and gas

         M_VAP = PC%ADJUST_EVAPORATION*M_VAP

         ! Add vapor or fuel gas to the grid cell

         IF (IGAS>0) YY(II,JJ,KK,IGAS)= MIN(1._EB,(WGT*M_VAP+M_GAS*Y_GAS)/(WGT*M_VAP+M_GAS))

         ! Add new mass from vaporized droplet to the grid cell

         RHO(II,JJ,KK) = RHO(II,JJ,KK) + WGT*M_VAP*RVC

         ! Get out of the loop if the droplet has evaporated completely

         IF (DR%R<=0._EB) CYCLE DROPLET_LOOP

         ! Assign liquid mass to the cell for airborne drops

         IF (DR%IOR==0) THEN
            DEN_ADD = WGT*M_DROP*RVC
            DROP_DEN(II,JJ,KK) = DROP_DEN(II,JJ,KK) + DEN_ADD
            DROP_TMP(II,JJ,KK) = DROP_TMP(II,JJ,KK) + DEN_ADD*DR%TMP
            DROP_RAD(II,JJ,KK) = DROP_RAD(II,JJ,KK) + DEN_ADD*DR%R
         ENDIF
         
      ENDDO DROPLET_LOOP

      ! Calculate cumulative quantities for droplet "clouds"

      WGT   = .5_EB
      OMWGT = 1._EB-WGT
    
      DROP_RAD = DROP_RAD/(DROP_DEN+1.E-10_EB)
      DROP_TMP = DROP_TMP/(DROP_DEN+1.E-10_EB)
    
      AVG_DROP_RAD(:,:,:,EVAP_INDEX ) = (AVG_DROP_DEN(:,:,:,EVAP_INDEX )*AVG_DROP_RAD(:,:,:,EVAP_INDEX )+DROP_DEN*DROP_RAD) &
                                          /(AVG_DROP_DEN(:,:,:,EVAP_INDEX ) + DROP_DEN + 1.0E-10_EB)
      AVG_DROP_TMP(:,:,:,EVAP_INDEX ) = (AVG_DROP_DEN(:,:,:,EVAP_INDEX )*AVG_DROP_TMP(:,:,:,EVAP_INDEX )+DROP_DEN*DROP_TMP) &
                                          /(AVG_DROP_DEN(:,:,:,EVAP_INDEX ) + DROP_DEN + 1.0E-10_EB)
      AVG_DROP_TMP(:,:,:,EVAP_INDEX ) = MAX(TMPM,AVG_DROP_TMP(:,:,:,EVAP_INDEX ))
      AVG_DROP_DEN(:,:,:,EVAP_INDEX ) = WGT*AVG_DROP_DEN(:,:,:,EVAP_INDEX ) + OMWGT*DROP_DEN
      WHERE (AVG_DROP_DEN(:,:,:,EVAP_INDEX )<0.0001_EB .AND. DROP_DEN==0._EB) AVG_DROP_DEN(:,:,:,EVAP_INDEX ) = 0.0_EB

   ENDDO PART_CLASS_LOOP

ENDDO EVAP_INDEX_LOOP

! Remove droplets that have completely evaporated
 
CALL REMOVE_DROPLETS(T,NM)
 
END SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER



SUBROUTINE PARTICLE_MOMENTUM_TRANSFER

! Add droplet momentum as a force term in momentum equation

USE COMP_FUNCTIONS, ONLY : SECOND
REAL(EB), POINTER, DIMENSION(:,:,:) :: FVXS,FVYS,FVZS
REAL(EB) :: FLUXMIN,FLUXMAX,XI,YJ,ZK
INTEGER :: II,JJ,KK,IIX,JJY,KKZ,I,J,K,IC,IW


FVXS  => WORK1
FVYS  => WORK2
FVZS  => WORK3
FVXS  = 0._EB
FVYS  = 0._EB
FVZS  = 0._EB
 
SUM_MOMENTUM_LOOP: DO I=1,NLP
   DR=>DROPLET(I)
   PC=>PARTICLE_CLASS(DR%CLASS)
   IF (DR%IOR/=0)   CYCLE SUM_MOMENTUM_LOOP
   IF (PC%MASSLESS) CYCLE SUM_MOMENTUM_LOOP
   XI  = CELLSI(FLOOR((DR%X-XS)*RDXINT))
   YJ  = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
   ZK  = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
   II  = FLOOR(XI+1._EB)
   JJ  = FLOOR(YJ+1._EB)
   KK  = FLOOR(ZK+1._EB)
   IF (SOLID(CELL_INDEX(II,JJ,KK))) CYCLE SUM_MOMENTUM_LOOP
   IIX = FLOOR(XI+.5_EB)
   JJY = FLOOR(YJ+.5_EB)
   KKZ = FLOOR(ZK+.5_EB)
   IC = CELL_INDEX(IIX,JJ,KK)
   IW = WALL_INDEX(IC,1)
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) FVXS(IIX,JJ,KK) = FVXS(IIX,JJ,KK) - DR%A_X
   IC = CELL_INDEX(II,JJY,KK)
   IW = WALL_INDEX(IC,2) 
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) FVYS(II,JJY,KK) = FVYS(II,JJY,KK) - DR%A_Y
   IC = CELL_INDEX(II,JJ,KKZ)
   IW = WALL_INDEX(IC,3) 
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) FVZS(II,JJ,KKZ) = FVZS(II,JJ,KKZ) - DR%A_Z
ENDDO SUM_MOMENTUM_LOOP
 
FLUXMAX =  100._EB
FLUXMIN = -FLUXMAX
 
DO K=0,KBAR
   DO J=0,JBAR
      DO I=0,IBAR
         IF (FVXS(I,J,K)>FLUXMAX) FVXS(I,J,K) = FLUXMAX
         IF (FVXS(I,J,K)<FLUXMIN) FVXS(I,J,K) = FLUXMIN
         FVX(I,J,K) = FVX(I,J,K) + FVXS(I,J,K)
         IF (FVYS(I,J,K)>FLUXMAX) FVYS(I,J,K) = FLUXMAX
         IF (FVYS(I,J,K)<FLUXMIN) FVYS(I,J,K) = FLUXMIN
         FVY(I,J,K) = FVY(I,J,K) + FVYS(I,J,K)
         IF (FVZS(I,J,K)>FLUXMAX) FVZS(I,J,K) = FLUXMAX
         IF (FVZS(I,J,K)<FLUXMIN) FVZS(I,J,K) = FLUXMIN
         FVZ(I,J,K) = FVZ(I,J,K) + FVZS(I,J,K)
      ENDDO
   ENDDO
ENDDO
 
END SUBROUTINE PARTICLE_MOMENTUM_TRANSFER
 
 
 
SUBROUTINE REMOVE_DROPLETS(T,NM)
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_DROPLETS 
! Remove droplets that do not lie in any mesh
 
INTEGER :: IKILL,I,NM,IPC
REAL(EB) :: T
REAL(EB), PARAMETER :: RDMIN=1.0E-6_EB   ! Minimum droplet size (m)
 
IKILL = 0
DROP_LOOP: DO I=1,NLP
   WEED_LOOP: DO
      DR=>DROPLET(I) 
      IPC=DR%CLASS 
      PC=>PARTICLE_CLASS(IPC)
      IF (I>NLP-IKILL) EXIT DROP_LOOP
      IF (DR%R<=RDMIN .AND. .NOT.PC%MASSLESS) THEN 
         CALL REPLACE
         CYCLE WEED_LOOP
      ENDIF
      IF (T-DR%T>PC%LIFETIME) THEN
         CALL REPLACE
         CYCLE WEED_LOOP
      ENDIF
      IF (DR%X>XS .AND. DR%X<XF .AND.DR%Y>YS .AND. DR%Y<YF .AND. DR%Z>ZS .AND. DR%Z<ZF)  CYCLE DROP_LOOP
      CALL REPLACE
   ENDDO WEED_LOOP
ENDDO DROP_LOOP
 
NLP = NLP - IKILL
 
CONTAINS
 
SUBROUTINE REPLACE
 
INTEGER :: OM,NOM
TYPE (MESH_TYPE), POINTER :: M
TYPE (OMESH_TYPE), POINTER :: M2
 
NOM = 0
SEARCH_LOOP: DO OM=1,NMESHES
   IF (NIC(NM,OM)==0) CYCLE SEARCH_LOOP
   IF (EVACUATION_ONLY(OM)) CYCLE SEARCH_LOOP
   M=>MESHES(OM)
   IF (DR%X>M%XS .AND. DR%X<M%XF .AND.  DR%Y>M%YS .AND. DR%Y<M%YF .AND.  DR%Z>M%ZS .AND. DR%Z<M%ZF) THEN
      NOM = OM
      EXIT SEARCH_LOOP
   ENDIF
ENDDO SEARCH_LOOP
 
IF (NOM/=0) THEN
   M2=>MESHES(NM)%OMESH(NOM)
   M2%N_DROP_ORPHANS = M2%N_DROP_ORPHANS + 1
   IF (M2%N_DROP_ORPHANS>M2%N_DROP_ORPHANS_DIM) CALL RE_ALLOCATE_DROPLETS(2,NM,NOM,1000)
   M2%DROPLET(M2%N_DROP_ORPHANS) = DROPLET(I)
ENDIF
 
DROPLET(I) = DROPLET(NLP-IKILL)
IKILL = IKILL + 1
 
END SUBROUTINE REPLACE
 
END SUBROUTINE REMOVE_DROPLETS
 
SUBROUTINE GET_REV_part(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') partrev(INDEX(partrev,':')+1:LEN_TRIM(partrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') partdate

END SUBROUTINE GET_REV_part
END MODULE PART