part.f90

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

MODULE PART
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE TRAN
USE TYPES, ONLY: DROPLET_TYPE, PARTICLE_CLASS_TYPE, PARTICLE_CLASS, WALL_TYPE
IMPLICIT NONE
PRIVATE
PUBLIC INSERT_DROPLETS_AND_PARTICLES,UPDATE_PARTICLES, &
       INITIALIZE_DROPLETS,INITIALIZE_TREES,GET_REV_part
TYPE (DROPLET_TYPE), POINTER :: DR
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC
CHARACTER(255), PARAMETER :: partid='$Id: part.f90 703 2007-09-27 22:33:46Z mcgratta $'
CHARACTER(255), PARAMETER :: partrev='$Revision: 703 $'
CHARACTER(255), PARAMETER :: partdate='$Date: 2007-09-27 18:33:46 -0400 (Thu, 27 Sep 2007) $'
 
CONTAINS
 

SUBROUTINE INITIALIZE_DROPLETS(NM)

! Insert droplets into the domain at the start of calculation

USE COMP_FUNCTIONS, ONLY : SECOND 
USE PHYSICAL_FUNCTIONS, ONLY : DROPLET_SIZE_DISTRIBUTION 
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_DROPLETS
 
REAL(EB) RN,MASS_SUM,LL,UL,BIN_SIZE,TNOW
REAL(EB) VOL1, VOL2, X1, X2, Y1, Y2, Z1, Z2
INTEGER I,J,II,JJ,KK,IL,IU,STRATUM,IPC
INTEGER, INTENT(IN) :: NM

IF (N_PART==0) RETURN            ! Don't waste time if no particles
IF (EVACUATION_ONLY(NM)) RETURN  ! Don't waste time if an evac mesh
 
TNOW=SECOND()
CALL POINT_TO_MESH(NM)
 
PART_CLASS_LOOP: DO IPC=1,N_PART
 
   PC=>PARTICLE_CLASS(IPC)

   ! Check if droplets evaporate into a SPECIES
 
   SPECIES_LOOP: DO I=1,N_SPECIES
      IF (PC%SPECIES=='null') CYCLE SPECIES_LOOP
      IF (SPECIES(I)%NAME==PC%SPECIES) THEN
         PC%SPECIES_INDEX = I
         EXIT SPECIES_LOOP
      ENDIF
   ENDDO SPECIES_LOOP
 
   ! If particles/droplets have a size distribution, initialize here
 
   IF_SIZE_DISTRIBUTION: IF (PC%DIAMETER > 0._EB) THEN
      CALL DROPLET_SIZE_DISTRIBUTION(PC%DIAMETER,PC%R_CDF(:),PC%CDF(:),NDC,PC%GAMMA,PC%SIGMA)
      BIN_SIZE = PC%R_CDF(NDC)/REAL(NSTRATA,EB)
      STRATIFY: DO I=1,NSTRATA
         LL = (I-1)*BIN_SIZE
         UL =  I   *BIN_SIZE
         LL_LOOP: DO J=1,NDC
            IF (PC%R_CDF(J)>LL) THEN
               IL = J-1 
               PC%IL_CDF(I) = J-1
               EXIT LL_LOOP
               ENDIF
            ENDDO LL_LOOP
         UL_LOOP: DO J=NDC,1,-1
            IF (PC%R_CDF(J)<=UL) THEN
               IU = J 
               PC%IU_CDF(I) = J
               EXIT UL_LOOP
               ENDIF
            ENDDO UL_LOOP
         PC%W_CDF(I) = PC%CDF(IU) - PC%CDF(IL)
      ENDDO STRATIFY
   ENDIF IF_SIZE_DISTRIBUTION
 
   ! If there is a specified number of initial droplets/particles, initialize
 
   IF (PC%TREE)         CYCLE PART_CLASS_LOOP
   IF (PC%N_INITIAL==0) CYCLE PART_CLASS_LOOP
 
   IF (PC%X1 == 0.0_EB .AND. PC%X2 == 0.0_EB .AND. PC%Y1 == 0.0_EB .AND. PC%Y2 == 0.0_EB .AND.  &
       PC%Z1 == 0.0_EB .AND. PC%Z2 == 0.0_EB) THEN
      X1 = XS 
      X2 = XF
      Y1 = YS 
      Y2 = YF
      Z1 = ZS 
      Z2 = ZF
      VOL2 = (XF - XS) * (YF - YS) * (ZF - ZS)
      VOL1 = VOL2
   ELSE
      IF (PC%X1>XF .OR. PC%X2<XS .OR. PC%Y1>YF .OR. PC%Y2<YS .OR. PC%Z1>ZF .OR. PC%Z2<ZS) CYCLE PART_CLASS_LOOP
      X1 = MAX(PC%X1,XS) 
      X2 = MIN(PC%X2,XF)
      Y1 = MAX(PC%Y1,YS) 
      Y2 = MIN(PC%Y2,YF)
      Z1 = MAX(PC%Z1,ZS) 
      Z2 = MIN(PC%Z2,ZF)
      VOL2 = (PC%X2 - PC%X1) * (PC%Y2 - PC%Y1) * (PC%Z2 - PC%Z1)
      VOL1 = (X2 - X1) * (Y2 - Y1) * (Z2 - Z1)
   ENDIF

   ! Assign properties to the initial droplets/particles
      
   MASS_SUM = 0._EB
   INITIALIZATION_LOOP: DO I=1,PC%N_INITIAL
      NLP = NLP + 1
      PARTICLE_TAG = PARTICLE_TAG + NMESHES
      IF (NLP>NLPDIM) THEN
         CALL RE_ALLOCATE_DROPLETS(1,NM,0,1000)
         DROPLET=>MESHES(NM)%DROPLET
      ENDIF
      DR=>DROPLET(NLP)
      BLOCK_OUT_LOOP:  DO
         CALL RANDOM_NUMBER(RN)
         DR%X = X1 + RN*(X2-X1)
         CALL RANDOM_NUMBER(RN)
         DR%Y = Y1 + RN*(Y2-Y1)
         CALL RANDOM_NUMBER(RN)
         DR%Z = Z1 + RN*(Z2-Z1)
         II = CELLSI(FLOOR((DR%X-XS)*RDXINT)) + 1._EB
         JJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT)) + 1._EB
         KK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT)) + 1._EB
         IF (.NOT.SOLID(CELL_INDEX(II,JJ,KK))) EXIT BLOCK_OUT_LOOP
      ENDDO BLOCK_OUT_LOOP
      DR%U   = 0._EB                     ! No initial velocity
      DR%V   = 0._EB
      DR%W   = 0._EB
      DR%TMP = PC%TMP_INITIAL            ! Initial temperature
      DR%T   = 0._EB                     ! Insertion time is 0
      DR%R   = 0._EB                     ! Radius is zero unless DIAMETER has been specified
      DR%PWT = 1._EB                     ! Weighting factor is one unless changed below
      DR%IOR = 0                         ! Orientation of solid surface (0 means the droplet/particle is not attached)
      DR%CLASS = IPC                     ! Class identifier
      DR%TAG   = PARTICLE_TAG            ! Unique integer tag
      IF (MOD(NLP,PC%SAMPLING)==0) THEN  ! Decide whether to show or not show in Smokeview
         DR%SHOW = .TRUE.    
      ELSE
         DR%SHOW = .FALSE.    
      ENDIF
 
      IF (PC%DIAMETER>0._EB) THEN
         IF (PC%MONODISPERSE) THEN
            DR%R   = 0.5_EB*PC%DIAMETER
            DR%PWT = 1._EB
         ELSE
            STRATUM = MOD(NLP-1,NSTRATA) + 1
            IL = PC%IL_CDF(STRATUM)
            IU = PC%IU_CDF(STRATUM)
            CALL RANDOM_CHOICE(PC%CDF(IL:IU),PC%R_CDF(IL:IU),IU-IL,DR%R)
            DR%PWT = PC%W_CDF(STRATUM)
            IF (2._EB*DR%R > PC%MAXIMUM_DIAMETER) THEN
               DR%PWT = DR%PWT*DR%R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
               DR%R = 0.5_EB*PC%MAXIMUM_DIAMETER
            ENDIF
         ENDIF
         MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3
      ENDIF
 
   ENDDO INITIALIZATION_LOOP
 
   ! Adjust particle weighting factor PWT so that desired MASS_PER_VOLUME is achieved

   IF (PC%DIAMETER>0._EB) DROPLET(NLP-PC%N_INITIAL+1:NLP)%PWT = &
      DROPLET(NLP-PC%N_INITIAL+1:NLP)%PWT*PC%MASS_PER_VOLUME*VOL1/MASS_SUM

ENDDO PART_CLASS_LOOP

TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE INITIALIZE_DROPLETS
 
 

SUBROUTINE INITIALIZE_TREES(NM)

USE MEMORY_FUNCTIONS, ONLY: RE_ALLOCATE_DROPLETS 
REAL(EB) DELTAZ_BIN,CANOPY_LENGTH,CANOPY_VOLUME,COSINE,TANGENT, &
         XC_MAX_ZC_MAX,XC_MAX_ZC_MIN,XC_MIN,XC_MAX, &
         YC_MIN,YC_MAX,ZBIN_VOLUME,ZC_MIN,ZC_MAX,CANOPY_WIDTH
INTEGER NCT,NNLP,NLP_TREE,NZB,NZBINS,P_PER_ZBIN,IPC
REAL(EB) RN,MASS_SUM
INTEGER I,II,JJ,KK,IL,IU,STRATUM
INTEGER, INTENT(IN) :: NM
 
CALL POINT_TO_MESH(NM)
 
TREE_LOOP: DO NCT=1,N_TREES
 
   IF (TREE_MESH(NCT)/=NM) CYCLE TREE_LOOP
   IPC = TREE_PARTICLE_CLASS(NCT)
   PC=>PARTICLE_CLASS(IPC)
 
   ! Build a conical tree
 
   CANOPY_WIDTH  = CANOPY_W(NCT)
   CANOPY_LENGTH = TREE_H(NCT) - CANOPY_B_H(NCT)
   TANGENT = 0.5_EB*CANOPY_W(NCT)/CANOPY_LENGTH
   COSINE = CANOPY_LENGTH/SQRT(CANOPY_LENGTH**2 + (0.5_EB*CANOPY_WIDTH)**2)
   NZBINS = 1._EB/(1._EB-COSINE)
   DELTAZ_BIN = CANOPY_LENGTH/REAL(NZBINS,EB)
   CANOPY_VOLUME = PI*CANOPY_WIDTH**2*CANOPY_LENGTH/12.
 
   NLP_TREE = 0
   ZBIN_LOOP: DO NZB=1,NZBINS
      ZC_MIN = CANOPY_B_H(NCT) + (NZB-1)*DELTAZ_BIN
      ZC_MAX = ZC_MIN + DELTAZ_BIN
      XC_MAX_ZC_MIN = (CANOPY_LENGTH-(ZC_MIN-CANOPY_B_H(NCT)))*TANGENT
      XC_MAX_ZC_MAX = (CANOPY_LENGTH-(ZC_MAX-CANOPY_B_H(NCT)))*TANGENT
      ZBIN_VOLUME = PI*DELTAZ_BIN/3._EB*(XC_MAX_ZC_MIN**2 + XC_MAX_ZC_MAX**2 +XC_MAX_ZC_MIN*XC_MAX_ZC_MAX )
      P_PER_ZBIN  = PC%N_INITIAL*ZBIN_VOLUME/CANOPY_VOLUME
 
      NISP_LOOP1: DO NNLP=1,P_PER_ZBIN
         NLP  = NLP + 1
         PARTICLE_TAG = PARTICLE_TAG + NMESHES
         NLP_TREE = NLP_TREE + 1
         IF (NLP>NLPDIM) THEN
            CALL RE_ALLOCATE_DROPLETS(1,NM,0,1000)
            DROPLET=>MESHES(NM)%DROPLET
            ENDIF
         DR=>DROPLET(NLP)
         DR%TAG = PARTICLE_TAG
         BLK_LOOP:  DO
            CALL RANDOM_NUMBER(RN)
            DR%Z = ZC_MIN + RN*(ZC_MAX-ZC_MIN)
            CALL RANDOM_NUMBER(RN)
            XC_MAX = (CANOPY_LENGTH-(DR%Z-CANOPY_B_H(NCT)))*TANGENT
            XC_MIN = -XC_MAX
            DR%X = XC_MIN + RN*(XC_MAX-XC_MIN)
            CALL RANDOM_NUMBER(RN)
            YC_MAX = SQRT(XC_MAX**2 - DR%X**2)
            YC_MIN = -YC_MAX
            DR%Y = YC_MIN + RN*(YC_MAX-YC_MIN)
            DR%X = X_TREE(NCT) + DR%X
            DR%Y = Y_TREE(NCT) + DR%Y
            DR%Z = Z_TREE(NCT) + DR%Z
            II = CELLSI(FLOOR((DR%X - XS)*RDXINT)) + 1._EB
            JJ = CELLSJ(FLOOR((DR%Y - YS)*RDYINT)) + 1._EB
            KK = CELLSK(FLOOR((DR%Z - ZS)*RDZINT)) + 1._EB
            IF (.NOT.SOLID(CELL_INDEX(II,JJ,KK))) EXIT BLK_LOOP
         ENDDO BLK_LOOP
      ENDDO NISP_LOOP1
   ENDDO ZBIN_LOOP
 
! Define physical characteristics of fuel elements
 
   MASS_SUM = 0._EB
   NISP_LOOP2: DO I=NLP-NLP_TREE+1,NLP
      DR=>DROPLET(I)
      DR%U = 0._EB
      DR%V = 0._EB
      DR%W = 0._EB
      DR%TMP = PC%TMP_INITIAL
      DR%T   = 0._EB
      DR%IOR = 0
      DR%A_X = 0._EB
      DR%A_Y = 0._EB
      DR%A_Z = 0._EB
      DR%CLASS = IPC
      IF (MOD(I,PC%SAMPLING)==0) THEN
         DR%SHOW = .TRUE.
      ELSE
         DR%SHOW = .FALSE.
      ENDIF
      IF (PC%MONODISPERSE) THEN
         DR%R   = 0.5_EB*PC%DIAMETER
         DR%PWT = 1._EB
      ELSE
         STRATUM = MOD(I-1,NSTRATA) + 1
         IL = PC%IL_CDF(STRATUM)
         IU = PC%IU_CDF(STRATUM)
         CALL RANDOM_CHOICE(PC%CDF(IL:IU),PC%R_CDF(IL:IU),IU-IL,DR%R)
         DR%PWT = PC%W_CDF(STRATUM)
         IF (2._EB*DR%R > PC%MAXIMUM_DIAMETER) THEN
            DR%PWT = DR%PWT*DR%R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
            DR%R = 0.5_EB*PC%MAXIMUM_DIAMETER
         ENDIF
      ENDIF
      MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3
   ENDDO NISP_LOOP2
 
   ! Adjust particle weighting factor PWT so that desired MASS_PER_VOLUME is achieved
 
   DROPLET(NLP-NLP_TREE+1:NLP)%PWT = DROPLET(NLP-NLP_TREE+1:NLP)%PWT* PC%MASS_PER_VOLUME*CANOPY_VOLUME/MASS_SUM
 
ENDDO TREE_LOOP
 
END SUBROUTINE INITIALIZE_TREES
 
 

SUBROUTINE INSERT_DROPLETS_AND_PARTICLES(T,NM)

! Insert sprinkler droplets and lagrangian particles into the domain every DT_INSERT seconds

USE COMP_FUNCTIONS, ONLY : SECOND 
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_DROPLETS
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE DEVICE_VARIABLES
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
REAL(EB) :: PHI_RN,RN,FLOW_RATE,THETA_RN,SPHI,CPHI,MASS_SUM, &
            STHETA,CTHETA,PWT0,DROPLET_SPEED,XI,YJ,ZK,SHIFT1,SHIFT2,XTMP,YTMP,ZTMP,VLEN, &
            TRIGT1,TRIGT2,TNOW,SOLID_ANGLE_FLOWRATE,SOLID_ANGLE_TOTAL_FLOWRATE,TSI
REAL(EB), PARAMETER :: VENT_OFFSET=0.1
INTEGER :: I,KS,II,JJ,KK,IC,IL,IU,STRATUM,IPC,DROP_SUM,IIG,JJG,KKG,IW,IBC,IOR
LOGICAL :: INSERT_ANOTHER_BATCH
TYPE (PROPERTY_TYPE), POINTER :: PY
TYPE (TABLES_TYPE), POINTER :: TA
TYPE (DEVICE_TYPE), POINTER :: DV
TYPE (SURFACE_TYPE), POINTER :: SF
 
IF (N_PART==0) RETURN  ! Don't waste time if no particles
TNOW=SECOND()
CALL POINT_TO_MESH(NM)

OVERALL_INSERT_LOOP: DO  ! Add more than one batch of particles/droplets if FDS time step is large enough
INSERT_ANOTHER_BATCH = .FALSE.

SPRINKLER_INSERT_LOOP: DO KS=1,N_DEVC  ! Loop over all devices, but look for sprinklers or nozzles
   DV => DEVICE(KS)
   PY => PROPERTY(DV%PROP_INDEX)
   IF (PY%PART_ID == 'null')   CYCLE SPRINKLER_INSERT_LOOP
   IF (DV%MESH/=NM)            CYCLE SPRINKLER_INSERT_LOOP
   IF (.NOT. DV%CURRENT_STATE) CYCLE SPRINKLER_INSERT_LOOP
   IPC = PROPERTY(DV%PROP_INDEX)%PART_INDEX
   PC=>PARTICLE_CLASS(IPC)
   IF (DV%T_CHANGE == T) THEN
      DV%T = T
      CYCLE SPRINKLER_INSERT_LOOP
   ENDIF
   IF (T < PC%INSERT_CLOCK(NM)) CYCLE SPRINKLER_INSERT_LOOP
   IF (T_BEGIN == DV%T_CHANGE) THEN
      TSI = T
   ELSE
      TSI = T - DV%T_CHANGE
   ENDIF
   FLOW_RATE = EVALUATE_RAMP(TSI,PY%FLOW_TAU,PY%FLOW_RAMP_INDEX)*PY%FLOW_RATE*(PC%DENSITY/1000._EB)/60._EB  ! kg/s
   IF (FLOW_RATE <= 0._EB) THEN
      DV%T = T
      CYCLE SPRINKLER_INSERT_LOOP
   ENDIF

   ! Direction initialization stuff

   TRIGT1 = ACOS(-DV%ORIENTATION(3))
   IF (DV%ORIENTATION(2)==0._EB) THEN
      TRIGT2 = ACOS(1._EB)
   ELSE
      TRIGT2 = ACOS(ABS(DV%ORIENTATION(1))/SQRT(DV%ORIENTATION(1)**2+DV%ORIENTATION(2)**2))
   ENDIF
 
   ! Droplet insertion loop
 
   MASS_SUM = 0._EB
   DROP_SUM = 0
   SOLID_ANGLE_TOTAL_FLOWRATE = 0._EB

   DROPLET_INSERT_LOOP: DO I=1,PC%N_INSERT

      IF (NLP+1>MAXIMUM_DROPLETS) EXIT DROPLET_INSERT_LOOP
      NLP = NLP+1
      PARTICLE_TAG = PARTICLE_TAG + NMESHES
      IF (NLP>NLPDIM) THEN
         CALL RE_ALLOCATE_DROPLETS(1,NM,0,1000)
         DROPLET=>MESHES(NM)%DROPLET
      ENDIF
      DR=>DROPLET(NLP)
 
      ! Set droplet propeties
 
      DR%TMP    = PC%TMP_INITIAL         ! Initial temperature
      DR%T      = T                      ! Time of insertion
      DR%IOR    = 0                      ! Orientation index (0 means it is not attached to a solid)
      DR%A_X    = 0._EB                  ! x component of drag on the gas in units of m/s2
      DR%A_Y    = 0._EB                  ! y component of drag on the gas in units of m/s2
      DR%A_Z    = 0._EB                  ! z component of drag on the gas in units of m/s2
      DR%CLASS = IPC                     ! The class of particles
      DR%TAG    = PARTICLE_TAG           ! A unique identifying integer
      IF (MOD(NLP,PC%SAMPLING)==0) THEN  ! Decide whether to show the droplet in Smokeview
         DR%SHOW = .TRUE.    
      ELSE
         DR%SHOW = .FALSE.
      ENDIF
 
      ! Randomly choose theta and phi
 
      CHOOSE_COORDS: DO
         PICK_PATTERN: IF(PY%SPRAY_PATTERN_INDEX>0) THEN !Use spray pattern table
            TA => TABLES(PY%SPRAY_PATTERN_INDEX)
            CALL RANDOM_NUMBER(RN)
            FIND_ROW: DO II=1,TA%NUMBER_ROWS
               IF (RN>PY%TABLE_ROW(II)) CYCLE FIND_ROW
               EXIT FIND_ROW
            END DO FIND_ROW
            CALL RANDOM_NUMBER(RN)
            THETA_RN = TA%TABLE_DATA(II,1) + RN*(TA%TABLE_DATA(II,2)-TA%TABLE_DATA(II,1))
            CALL RANDOM_NUMBER(RN)
            PHI_RN = TA%TABLE_DATA(II,3) + RN*(TA%TABLE_DATA(II,4)-TA%TABLE_DATA(II,3))
            SOLID_ANGLE_FLOWRATE = TA%TABLE_DATA(II,6)
            SOLID_ANGLE_TOTAL_FLOWRATE = SOLID_ANGLE_TOTAL_FLOWRATE + TA%TABLE_DATA(II,6)
            PY%DROPLET_VELOCITY  = TA%TABLE_DATA(II,5)
         ELSE PICK_PATTERN !Use conical spray