part.f90

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

   VREL  = DR%V - VBAR
   WREL  = DR%W - WBAR
   QREL  = MAX(1.E-6_EB,SQRT(UREL*UREL + VREL*VREL + WREL*WREL))
   TMP_G = MAX(TMPMIN,TMP(II,JJ,KK))
   RHO_G = RHO(II,JJ,KK)
   MU_AIR = SPECIES(0)%MU(MIN(500,NINT(0.1_EB*TMP_G)))
   DR%RE  = RHO_G*QREL*2._EB*RD/MU_AIR
   DRAG_CALC: IF (DR%IOR==0 .AND. DR%RE>0) THEN
      C_DRAG  = DRAG(DR%RE)
      SFAC    = DR%PWT*RDS*PIO2*QREL*C_DRAG
      DR%A_X  = SFAC*UREL*RVC
      DR%A_Y  = SFAC*VREL*RVC
      DR%A_Z  = SFAC*WREL*RVC
      IF (.NOT.PC%STATIC) THEN
         CONST   = 8._EB*PC%DENSITY*RD/(3._EB*RHO_G*C_DRAG*QREL)
         BFAC    = EXP(-DT/CONST)
         AUREL   = CONST*GRVT1
         AVREL   = CONST*GRVT2
         AWREL   = CONST*GRVT3
         DR%U    = UBAR + (UREL+AUREL)*BFAC - AUREL
         DR%V    = VBAR + (VREL+AVREL)*BFAC - AVREL
         DR%W    = WBAR + (WREL+AWREL)*BFAC - AWREL
      ENDIF
   ELSE
      DR%A_X  = 0._EB
      DR%A_Y  = 0._EB
      DR%A_Z  = 0._EB
   ENDIF DRAG_CALC

   ! If the particle does not move, but does drag, go on to the next particle

   IF (PC%STATIC) CYCLE DROPLET_LOOP
 
   ! Decide how many time steps to use in tracking particle
    
   DTMIN = DT
   UVW = MAX( ABS(DR%U)*RDX(II),ABS(DR%V)*RDY(JJ),ABS(DR%W)*RDZ(KK) )
   IF (UVW/=0._EB) DTMIN = MIN(DTMIN,1._EB/UVW)
    
   NITER = 1
   DTSP  = DT
   NLOOP: DO N=0,3
      IF (DTMIN<DT*0.5_EB**N) THEN
         NITER = 2**(N+1)
         DTSP =  DT*0.5_EB**(N+1)
      ENDIF
   ENDDO NLOOP
    
   ! Iterate over a single time step
    
   SUB_TIME_STEP_ITERATIONS: DO N=1,NITER
    
      ! Get current particle coordinates

      IF (N>1) THEN
         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)
         IC = CELL_INDEX(II,JJ,KK)
      ENDIF
    
      ! Update droplet position

      X_OLD = DR%X
      Y_OLD = DR%Y
      Z_OLD = DR%Z
      DR%X  = DR%X + DR%U*DTSP
      DR%Y  = DR%Y + DR%V*DTSP
      DR%Z  = DR%Z + DR%W*DTSP
    
      ! Droplet hits the floor
    
      IF (POROUS_FLOOR .AND. DR%Z<ZS .AND. PC%EVAP_INDEX>0) THEN
         IC = CELL_INDEX(II,JJ,1)
         IW = WALL_INDEX(IC,-3)
         IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY .AND. ACCUMULATE_WATER) THEN
            AWMPUA(IW,PC%EVAP_INDEX) = AWMPUA(IW,PC%EVAP_INDEX) + DR%PWT*PC%FTPR*RDC*RAW(IW)
         ENDIF
         CYCLE DROPLET_LOOP
      ENDIF
    
      IF (.NOT.POROUS_FLOOR .AND. DR%Z<ZS) THEN
         IC = CELL_INDEX(II,JJ,1)
         IW = WALL_INDEX(IC,-3)
         IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
            IF (ACCUMULATE_WATER) AWMPUA(IW,PC%EVAP_INDEX) = AWMPUA(IW,PC%EVAP_INDEX) + DR%PWT*PC%FTPR*RDC*RAW(IW)
            DR%Z = ZS + 0.05_EB*DZ(1)
            DR%IOR = 3
            CALL RANDOM_NUMBER(RN)
            THETA_RN = TWOPI*RN
            DR%U = PC%HORIZONTAL_VELOCITY*COS(THETA_RN)
            DR%V = PC%HORIZONTAL_VELOCITY*SIN(THETA_RN)
            DR%W = 0._EB
            ENDIF
      ENDIF
    
      ! Where is the droplet now?
    
      XI  = CELLSI(FLOOR((DR%X-XS)*RDXINT))
      YJ  = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
      ZK  = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
      IIN = FLOOR(XI+1._EB)
      JJN = FLOOR(YJ+1._EB)
      KKN = FLOOR(ZK+1._EB)
      IF (IIN<0 .OR. IIN>IBP1) CYCLE DROPLET_LOOP
      IF (JJN<0 .OR. JJN>JBP1) CYCLE DROPLET_LOOP
      IF (KKN<0 .OR. KKN>KBP1) CYCLE DROPLET_LOOP
      ICN = CELL_INDEX(IIN,JJN,KKN)

      IF (IC==0 .OR. ICN==0) CYCLE SUB_TIME_STEP_ITERATIONS

      IF (.NOT.SOLID(ICN)) THEN
         IF (DR%X<XS .OR. DR%X>XF) CYCLE DROPLET_LOOP
         IF (DR%Y<YS .OR. DR%Y>YF) CYCLE DROPLET_LOOP
         IF (DR%Z<ZS .OR. DR%Z>ZF) CYCLE DROPLET_LOOP
      ENDIF

      ! If droplet hits an obstacle, change its properties

      AIR_TO_SOLID: IF (II/=IIN .OR. JJ/=JJN .OR. KK/=KKN) THEN

         IOR_OLD   = DR%IOR
         HIT_SOLID = .FALSE.

         ! Check if any solid boundaries of original grid cell have been crossed

         IWP1 = WALL_INDEX(IC, 1) 
         IWM1 = WALL_INDEX(IC,-1)
         IWP2 = WALL_INDEX(IC, 2)
         IWM2 = WALL_INDEX(IC,-2)
         IWP3 = WALL_INDEX(IC, 3)
         IWM3 = WALL_INDEX(IC,-3)
         STEP_FRACTION = 1._EB

         IF (KKN>KK .AND. BOUNDARY_TYPE(IWP3)==SOLID_BOUNDARY) THEN
            DR%IOR=-3
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KK)-Z_OLD-0.05_EB*DZ(KK))/(DR%Z-Z_OLD))
         ENDIF
         IF (KKN<KK .AND. BOUNDARY_TYPE(IWM3)==SOLID_BOUNDARY) THEN
            DR%IOR= 3
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KKN)-Z_OLD+0.05_EB*DZ(KKN))/(DR%Z-Z_OLD))
         ENDIF
         IF (IIN>II .AND. BOUNDARY_TYPE(IWP1)==SOLID_BOUNDARY) THEN
            DR%IOR=-1
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(II)-X_OLD-0.05_EB*DX(II))/(DR%X-X_OLD))
         ENDIF
         IF (IIN<II .AND. BOUNDARY_TYPE(IWM1)==SOLID_BOUNDARY) THEN
            DR%IOR= 1
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(IIN)-X_OLD+0.05_EB*DX(IIN))/(DR%X-X_OLD))
         ENDIF
         IF (JJN>JJ .AND. BOUNDARY_TYPE(IWP2)==SOLID_BOUNDARY) THEN
            DR%IOR=-2
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJ)-Y_OLD-0.05_EB*DY(JJ))/(DR%Y-Y_OLD))
         ENDIF
         IF (JJN<JJ .AND. BOUNDARY_TYPE(IWM2)==SOLID_BOUNDARY) THEN
            DR%IOR= 2
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJN)-Y_OLD+0.05_EB*DY(JJN))/(DR%Y-Y_OLD))
         ENDIF

         IML = MINLOC(STEP_FRACTION,DIM=1)
         IOR_FIRST = 0
         SELECT CASE(IML)
            CASE(1)
               IOR_FIRST = -3
            CASE(2)
               IOR_FIRST = -2
            CASE(3)
               IOR_FIRST = -1
            CASE(5)
               IOR_FIRST =  1
            CASE(6)
               IOR_FIRST =  2
            CASE(7)
               IOR_FIRST =  3
         END SELECT
         DR%WALL_INDEX = WALL_INDEX(IC,-IOR_FIRST)
         
         ! If no solids boundaries of original cell have been crossed, check boundaries of new grid cell
    
         IF (DR%WALL_INDEX==0) THEN
            IWP1 = WALL_INDEX(ICN, 1)
            IWM1 = WALL_INDEX(ICN,-1)
            IWP2 = WALL_INDEX(ICN, 2)
            IWM2 = WALL_INDEX(ICN,-2)
            IWP3 = WALL_INDEX(ICN, 3)
            IWM3 = WALL_INDEX(ICN,-3)
            HIT_SOLID = .FALSE.
            STEP_FRACTION = 1._EB
            IF (KKN>KK .AND. BOUNDARY_TYPE(IWM3)==SOLID_BOUNDARY) THEN
               DR%IOR=-3
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KK)-Z_OLD-0.05_EB*DZ(KK))/(DR%Z-Z_OLD))
            ENDIF
            IF (KKN<KK .AND. BOUNDARY_TYPE(IWP3)==SOLID_BOUNDARY) THEN
               DR%IOR= 3
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KKN)-Z_OLD+0.05_EB*DZ(KKN))/(DR%Z-Z_OLD))
            ENDIF
            IF (IIN>II .AND. BOUNDARY_TYPE(IWM1)==SOLID_BOUNDARY) THEN
               DR%IOR=-1
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(II)-X_OLD-0.05_EB*DX(II))/(DR%X-X_OLD))
            ENDIF
            IF (IIN<II .AND. BOUNDARY_TYPE(IWP1)==SOLID_BOUNDARY) THEN
               DR%IOR= 1
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(IIN)-X_OLD+0.05_EB*DX(IIN))/(DR%X-X_OLD))
            ENDIF
            IF (JJN>JJ .AND. BOUNDARY_TYPE(IWM2)==SOLID_BOUNDARY) THEN
               DR%IOR=-2
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJ)-Y_OLD-0.05_EB*DY(JJ))/(DR%Y-Y_OLD))
            ENDIF
            IF (JJN<JJ .AND. BOUNDARY_TYPE(IWP2)==SOLID_BOUNDARY) THEN
               DR%IOR= 2
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJN)-Y_OLD+0.05_EB*DY(JJN))/(DR%Y-Y_OLD))
            ENDIF

            IML = MINLOC(STEP_FRACTION,DIM=1)
            IOR_FIRST = 0
            SELECT CASE(IML)
               CASE(1)
                  IOR_FIRST = -3
               CASE(2)
                  IOR_FIRST = -2
               CASE(3)
                  IOR_FIRST = -1
               CASE(5)
                  IOR_FIRST =  1
               CASE(6)
                  IOR_FIRST =  2
               CASE(7)
                  IOR_FIRST =  3
            END SELECT
            DR%WALL_INDEX = WALL_INDEX(ICN,IOR_FIRST)
         ENDIF
   
         ! Check if droplet has crossed no solid planes or too many
  
         IF_HIT_SOLID: IF (HIT_SOLID) THEN
 
            IF (DR%WALL_INDEX==0) CYCLE SUB_TIME_STEP_ITERATIONS

            ! Adjust the size of the droplet and weighting factor 

            R_NEW  = MIN(0.5_EB*SURFACE_DROPLET_DIAMETER,(DR%PWT*DR%R**3)**ONTH)
            DR%PWT = DR%PWT*DR%R**3/R_NEW**3
            DR%R   = R_NEW

            ! Move particle to where it almost hits solid

            DR%X = X_OLD + MINVAL(STEP_FRACTION)*DTSP*DR%U
            DR%Y = Y_OLD + MINVAL(STEP_FRACTION)*DTSP*DR%V
            DR%Z = Z_OLD + MINVAL(STEP_FRACTION)*DTSP*DR%W

            XI  = CELLSI(FLOOR((DR%X-XS)*RDXINT))
            YJ  = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
            ZK  = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
            IIN = FLOOR(XI+1._EB)
            JJN = FLOOR(YJ+1._EB)
            KKN = FLOOR(ZK+1._EB)
            ICN = CELL_INDEX(IIN,JJN,KKN)

            IF (IOR_OLD==DR%IOR) CYCLE SUB_TIME_STEP_ITERATIONS

            ! Choose a direction for the droplets to move

            DIRECTION: SELECT CASE(ABS(DR%IOR))
               CASE (1:2) DIRECTION  
                  DR%U = 0._EB
                  DR%V = 0._EB
                  DR%W = -PC%VERTICAL_VELOCITY 
               CASE(3) DIRECTION
                  CALL RANDOM_NUMBER(RN)
                  THETA_RN = TWOPI*RN
                  DR%U = PC%HORIZONTAL_VELOCITY*COS(THETA_RN)
                  DR%V = PC%HORIZONTAL_VELOCITY*SIN(THETA_RN)
                  DR%W = 0._EB
            END SELECT DIRECTION

         ENDIF IF_HIT_SOLID
      ENDIF AIR_TO_SOLID 

      ! Check if droplets that were attached to a solid are still attached after the time update

      IW = 0
      SELECT CASE(DR%IOR)
         CASE( 1)
            IW = WALL_INDEX(ICN,-1)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%X = DR%X - 0.2_EB*DX(II)
               DR%W = -DR%W
            ENDIF
         CASE(-1)
            IW = WALL_INDEX(ICN, 1)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%X = DR%X + 0.2_EB*DX(II)
               DR%W = -DR%W
            ENDIF
         CASE( 2)
            IW = WALL_INDEX(ICN,-2)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%Y = DR%Y - 0.2_EB*DY(JJ)
               DR%W = -DR%W
            ENDIF
         CASE(-2)
            IW = WALL_INDEX(ICN, 2)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%Y = DR%Y + 0.2_EB*DY(JJ)
               DR%W = -DR%W
            ENDIF
         CASE( 3)
            IW = WALL_INDEX(ICN,-3)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN  ! Particle has reached the edge of a horizontal surface
               DR%U = -DR%U
               DR%V = -DR%V
               DR%Z =  DR%Z - 0.2_EB*DZ(KK)
            ENDIF
         CASE(-3)
            IW = WALL_INDEX(ICN, 3)
      END SELECT

      IF (DR%IOR/=0 .AND. BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
         DR%IOR = 0
         DR%WALL_INDEX = 0
      ELSE
         DR%WALL_INDEX = WALL_INDEX(ICN,-DR%IOR)
      ENDIF

   ENDDO SUB_TIME_STEP_ITERATIONS

ENDDO DROPLET_LOOP

! Remove out-of-bounds particles
 
CALL REMOVE_DROPLETS(T,NM)

END SUBROUTINE MOVE_PARTICLES



SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,NM)
    
! Mass and energy transfer between gas and droplets

USE PHYSICAL_FUNCTIONS, ONLY : GET_MASS_FRACTION2

REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_DEN,DROP_RAD,DROP_TMP
REAL(EB), POINTER, DIMENSION(:) :: FILM_THICKNESS
REAL(EB) :: R_DROP,NUSSELT,K_AIR,H_V, &
            RVC,WGT,OMWGT,Q_CON_GAS,Q_CON_WALL,Q_RAD,H_HEAT,H_MASS,SH_FAC_GAS,SH_FAC_WALL,NU_FAC_GAS,NU_FAC_WALL, &
            T,PR_AIR,M_VAP,XI,YJ,ZK,RDT,MU_AIR,H_SOLID,Q_DOT_RAD,DEN_ADD, &
            Y_DROP,Y_GAS,LENGTH,U2,V2,W2,VEL,DENOM,DY_DTMP_DROP,TMP_DROP_NEW,TMP_WALL,H_WALL, &
            SC_AIR,D_AIR,DHOR,SHERWOOD,X_DROP,M_DROP,RHO_G,MW_RATIO,MW_DROP,FTPR,&
            C_DROP,M_GAS,A_DROP,TMP_G,TMP_DROP,TMP_MELT,TMP_BOIL,MINIMUM_FILM_THICKNESS,RE_L,OMRAF,Q_TEMP
INTEGER :: I,II,JJ,KK,IW,IGAS,N_PC,EVAP_INDEX
INTEGER, INTENT(IN) :: NM
REAL(EB), PARAMETER :: RUN_AVG_FAC=0.5_EB

! Initializations
RDT    = 1._EB/DT
OMRAF  = 1._EB - RUN_AVG_FAC

! Rough estimates

MINIMUM_FILM_THICKNESS = 1.E-5_EB   ! Minimum thickness of liquid film on the surface (m)
H_SOLID                = 300._EB    ! Heat transfer coefficient from solid surface to drop (W/m2/K)

! Empirical coefficients

D_AIR                  = 2.6E-5_EB  ! Water Vapor - Air binary diffusion (m2/s at 25 C, Incropera & DeWitt, Table A.8) 
SC_AIR                 = 0.6_EB     ! NU_AIR/D_AIR (Incropera & DeWitt, Chap 7, External Flow)
PR_AIR                 = 0.7_EB     
SH_FAC_GAS             = 0.6_EB*SC_AIR**ONTH
NU_FAC_GAS             = 0.6_EB*PR_AIR**ONTH        
SH_FAC_WALL            = 0.037_EB*SC_AIR**ONTH
NU_FAC_WALL            = 0.037_EB*PR_AIR**ONTH        

! Working arrays

IF (N_EVAP_INDICIES>0) THEN
   D_VAP  = 0._EB
   WCPUA  = RUN_AVG_FAC*WCPUA
   WMPUA  = RUN_AVG_FAC*WMPUA
   DROP_DEN => WORK4
   DROP_RAD => WORK5
   DROP_TMP => WORK6
ENDIF

! Loop over all types of evaporative species

EVAP_INDEX_LOOP: DO EVAP_INDEX = 1,N_EVAP_INDICIES

   FILM_THICKNESS => WALL_WORK2
   FILM_THICKNESS =  0._EB