init.f90

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

ENDIF
IF (I_OBST>0 .AND. SURFACE(IBCX)%MASS_FLUX_TOTAL/=-999._EB) THEN
   OBX=>M%OBSTRUCTION(I_OBST)
   M%UW0(IW) = SURFACE(IBCX)%MASS_FLUX_TOTAL / RHOA * M%AREA_ADJUST(IW)
ENDIF
 
T_ACTIVATE = -1._EB
VENT_FOUND = .FALSE.
 
! Check if there is a vent embedded in the surface
 
VLOOP: DO N=1,M%N_VENT
 
   VT => M%VENTS(N)
   IF (I_OBST>0) THEN
      IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY)       CYCLE VLOOP
      IF (.NOT.M%OBSTRUCTION(I_OBST)%ALLOW_VENT) CYCLE VLOOP
   ENDIF
   IF (VT%IOR/=IOR) CYCLE VLOOP
   IF (IBCX==INTERPOLATED_SURF_INDEX) CYCLE VLOOP
   IF ((VT%BOUNDARY_TYPE==OPEN_BOUNDARY .OR. VT%BOUNDARY_TYPE==MIRROR_BOUNDARY) .AND.M%IJKW(9,IW)>0) CYCLE VLOOP
 
   IF (ABS(IOR)==1) THEN
      IF (IOR== 1 .AND. I/=VT%I1  ) CYCLE VLOOP
      IF (IOR==-1 .AND. I/=VT%I1+1) CYCLE VLOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VLOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VLOOP
   ENDIF
   IF (ABS(IOR)==2) THEN
      IF (IOR== 2 .AND. J/=VT%J1  ) CYCLE VLOOP
      IF (IOR==-2 .AND. J/=VT%J1+1) CYCLE VLOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VLOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VLOOP
   ENDIF
   IF (ABS(IOR)==3) THEN
      IF (IOR== 3 .AND. K/=VT%K1  ) CYCLE VLOOP
      IF (IOR==-3 .AND. K/=VT%K1+1) CYCLE VLOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VLOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VLOOP
   ENDIF
 
   ! Check if there are over-lapping VENTs
  
   IF (VENT_FOUND) THEN
      WRITE(MESSAGE,'(A,I2,A,3I3,A,I3)') 'WARNING: Two VENTs overlap in MESH ',NM,', Cell',I,J,K,'. FDS rejects VENT ',VT%ORDINAL
      WRITE(0,'(/A/)') MESSAGE
      EXIT VLOOP
   ENDIF

   VENT_FOUND = .TRUE.

   ! Reassign the SURF index to be that of the VENT

   M%VENT_INDEX(IW) = N
   IBCX = VT%IBC
   M%AREA_ADJUST(IW) = VT%INPUT_AREA/VT%FDS_AREA
   IF (M%AREA_ADJUST(IW)==0._EB) M%AREA_ADJUST(IW) = 1._EB
 
   ! Set the velocity at each surface cell
 
   IF (SURFACE(IBCX)%VEL/=-999._EB) THEN
      M%UW0(IW) = SURFACE(IBCX)%VEL
   ELSE
      M%UW0(IW) = 0._EB
   ENDIF
 
   IF (SURFACE(IBCX)%VOLUME_FLUX    /=-999._EB)   M%UW0(IW) = SURFACE(IBCX)%VOLUME_FLUX/VT%TOTAL_INPUT_AREA
   IF (SURFACE(IBCX)%MASS_FLUX_TOTAL/=-999._EB)   M%UW0(IW) = SURFACE(IBCX)%MASS_FLUX_TOTAL/RHOA*M%AREA_ADJUST(IW)   
 
   ! Special velocity profiles
 
   IF (SURFACE(IBCX)%PROFILE==PARABOLIC) THEN 
      SELECT CASE(ABS(IOR))
         CASE(1)
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            Z1 = M%Z(VT%K1)
            Z2 = M%Z(VT%K2)
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            PZ = 4._EB*(M%ZC(K)-Z1)*(Z2-M%ZC(K))/(Z2-Z1)**2
            M%UW0(IW) = M%UW0(IW)*PY*PZ
         CASE(2)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            Z1 = M%Z(VT%K1)
            Z2 = M%Z(VT%K2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            PZ = 4._EB*(M%ZC(K)-Z1)*(Z2-M%ZC(K))/(Z2-Z1)**2
            M%UW0(IW) = M%UW0(IW)*PX*PZ
         CASE(3)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            IF (CYLINDRICAL .AND. X1==0._EB) X1 = -X2
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            IF (CYLINDRICAL) THEN
               M%UW0(IW) = M%UW0(IW)*PX
               ELSE
               M%UW0(IW) = M%UW0(IW)*PX*PY
               ENDIF
      END SELECT
   ENDIF
 
   IF (SURFACE(IBCX)%PROFILE==ATMOSPHERIC) THEN
      Z1 = M%Z(VT%K1)
      M%UW0(IW) = M%UW0(IW)*((M%ZC(K)-Z1)/  SURFACE(IBCX)%Z0)**SURFACE(IBCX)%PLE
   ENDIF
 
   IF (SURFACE(IBCX)%PROFILE==ONED_PARABOLIC) THEN 
      SELECT CASE(ABS(IOR))
         CASE(1)
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            M%UW0(IW) = M%UW0(IW)*PY
         CASE(2)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            M%UW0(IW) = M%UW0(IW)*PX
         CASE(3)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            M%UW0(IW) = M%UW0(IW)*PX*PY
      END SELECT
   ENDIF
 
   ! Check if fire spreads radially
 
   IF (VT%X0>-900._EB) THEN
      DIST = SQRT((M%XW(IW)-VT%X0)**2 +(M%YW(IW)-VT%Y0)**2 +(M%ZW(IW)-VT%Z0)**2)
      T_ACTIVATE = T_BEGIN+DIST/VT%FIRE_SPREAD_RATE
   ENDIF

   ! Miscellaneous settings
 
   IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY) THEN
      M%BOUNDARY_TYPE(IW) = OPEN_BOUNDARY
      IF (M%PRESSURE_ZONE_WALL(IW)>0) THEN
         WRITE(MESSAGE,'(A,I2,A)') 'ERROR: Pressure ZONE ',M%PRESSURE_ZONE_WALL(IW),' cannot have an OPEN boundary'
         CALL SHUTDOWN(MESSAGE)
      ENDIF
   ENDIF

   IF (VT%BOUNDARY_TYPE==MIRROR_BOUNDARY) M%BOUNDARY_TYPE(IW) = MIRROR_BOUNDARY
   IF (.NOT.VT%ACTIVATED) T_ACTIVATE = 1000000._EB
   M%IJKW(5,IW) = IBCX
 
ENDDO VLOOP
 
! Set ignition time of each boundary cell
 
IF (T_ACTIVATE<0._EB) THEN
   M%TW(IW) = SURFACE(IBCX)%T_IGN
ELSE
   M%TW(IW) = T_ACTIVATE
ENDIF
 
! Initialize solid temperature and densities
 
SF => SURFACE(IBCX)
IF (SF%THERMALLY_THICK) THEN
   WC => M%WALL(IW)
   ALLOCATE(WC%TMP_S(0:SF%N_CELLS+1),STAT=IZERO)
   CALL ChkMemErr('INIT','TMP_S',IZERO)
   WC%TMP_S(:)   = SF%TMP_INNER
   ALLOCATE(WC%RHO_S(0:SF%N_CELLS+1,SF%N_MATL),STAT=IZERO)
   CALL ChkMemErr('INIT','RHO_S',IZERO)
   WC%RHO_S = 0._EB
   DO II=0,SF%N_CELLS+1
      IL = SF%LAYER_INDEX(II)
      DO NN=1,SF%N_LAYER_MATL(IL)
         DO N=1,SF%N_MATL
            IF (SF%LAYER_MATL_INDEX(IL,NN)==SF%MATL_INDEX(N)) &
               WC%RHO_S(II,N) = SF%LAYER_MATL_FRAC(IL,NN)*SF%LAYER_DENSITY(IL)
         ENDDO
      ENDDO
   ENDDO
   IF (SF%SHRINK) THEN
      ALLOCATE(WC%LAYER_THICKNESS(SF%N_LAYERS),STAT=IZERO)
      CALL ChkMemErr('INIT','LAYER_THICKNESS',IZERO)
      WC%LAYER_THICKNESS = SF%LAYER_THICKNESS(1:SF%N_LAYERS)
      ALLOCATE(WC%N_LAYER_CELLS(SF%N_LAYERS),STAT=IZERO)
      CALL ChkMemErr('INIT','N_LAYER_CELLS',IZERO)
      WC%N_LAYER_CELLS = SF%N_LAYER_CELLS
      ALLOCATE(WC%X_S(0:SF%N_CELLS),STAT=IZERO)
      CALL ChkMemErr('INIT','X_S',IZERO)
      WC%X_S = SF%X_S
   ENDIF
   WC%SHRINKING = .FALSE.
   WC%BURNAWAY  = .FALSE.
   M%TMP_F(IW) = SF%TMP_INNER
   M%TMP_B(IW) = SF%TMP_INNER
   M%TMP_W(IW) = SF%TMP_INNER
ENDIF

 
! Misc
 
M%E_WALL(IW) = SF%EMISSIVITY
M%NPPCW(IW) = SF%NPPC      ! Number of particles per cell

END SUBROUTINE INIT_WALL_CELL
 
 
SUBROUTINE OPEN_AND_CLOSE(T,NM)

! Check to see if a cell or OBSTruction is to be created or removed, or a VENT activated of deactivated

USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_STRINGS
USE CONTROL_VARIABLES, ONLY : CONTROL
USE DEVICE_VARIABLES, ONLY : DEVICE
REAL(EB) :: T
INTEGER  :: N,II,JJ,KK,IW,IOR
INTEGER, INTENT(IN) :: NM
LOGICAL :: CREATE_OBST,REMOVE_OBST,ACTIVATE_VENT,DEACTIVATE_VENT
CHARACTER(12) :: SV_LABEL
TYPE (VENTS_TYPE), POINTER :: VT
 
CALL POINT_TO_MESH(NM)
 
! Check to see if an obstacle is to be removed or created
OBST_LOOP: DO N=1,N_OBST
   OB=>OBSTRUCTION(N)
   IF (.NOT. OB%REMOVABLE) CYCLE OBST_LOOP
   CREATE_OBST = .FALSE.
   REMOVE_OBST = .FALSE.
   IF (OB%CONSUMABLE .AND. OB%MASS <= 0._EB) REMOVE_OBST = .TRUE.

   IF (OB%DEVC_INDEX > 0) THEN
      IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
      IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
         CREATE_OBST = .TRUE.
      ELSE
         REMOVE_OBST = .TRUE.      
      ENDIF
   ELSEIF (OB%CTRL_INDEX > 0) THEN
      IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
      IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
         CREATE_OBST = .TRUE.
      ELSE
         REMOVE_OBST = .TRUE.      
      ENDIF
   ENDIF

  SV_LABEL  = 'null'

  IF (CREATE_OBST .AND. OB%HIDDEN) THEN
     OB%HIDDEN = .FALSE.
     SV_LABEL  = 'SHOW_OBST'
     CALL CREATE_OR_REMOVE_OBST(NM,OB%I1+1,OB%I2,OB%J1+1,OB%J2,OB%K1+1,OB%K2,1,N)
  ENDIF

   IF (REMOVE_OBST .AND. (.NOT. OB%HIDDEN)) THEN
      OB%HIDDEN = .TRUE.
      SV_LABEL  = 'HIDE_OBST'
      CALL CREATE_OR_REMOVE_OBST(NM,OB%I1+1,OB%I2,OB%J1+1,OB%J2,OB%K1+1,OB%K2,0,N)
   ENDIF

   ! Write a message to the Smokeview .smv file that the obstruction has been created or removed

   IF (SV_LABEL /= 'null') THEN
      IF (N_STRINGS+2>N_STRINGS_MAX) THEN
         CALL RE_ALLOCATE_STRINGS(NM)
         STRING => MESHES(NM)%STRING
      ENDIF
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(I6,F10.2)') N,T
   ENDIF

ENDDO OBST_LOOP
 
! Check to see if a vent should be activated or deactivated
 
VENT_LOOP: DO N=1,N_VENT
   VT => VENTS(N)

   ACTIVATE_VENT   = .FALSE.
   DEACTIVATE_VENT = .FALSE.

   IF (VT%DEVC_INDEX > 0) THEN
      IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(VT%DEVC_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
      IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE) THEN
         ACTIVATE_VENT   = .TRUE.
      ELSE
         DEACTIVATE_VENT = .TRUE.
      ENDIF
   ELSEIF (VT%CTRL_INDEX > 0) THEN
      IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(VT%CTRL_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
      IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE) THEN
         ACTIVATE_VENT   = .TRUE.
      ELSE
         DEACTIVATE_VENT = .TRUE.
      ENDIF
   ENDIF

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) CYCLE VENT_LOOP

   ! Look through all boundary cells looking for those that conform to the given VENT

   SEARCH: DO IW=1,NWC
      II  = IJKW(1,IW)
      JJ  = IJKW(2,IW)
      KK  = IJKW(3,IW)
      IOR = IJKW(4,IW)
      IF (IOR/=VT%IOR) CYCLE SEARCH
      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (IOR== 1 .AND.  II/=VT%I1)   CYCLE SEARCH
            IF (IOR==-1 .AND.  II/=VT%I1+1) CYCLE SEARCH
            IF (JJ<VT%J1+1 .OR. JJ>VT%J2)   CYCLE SEARCH
            IF (KK<VT%K1+1 .OR. KK>VT%K2)   CYCLE SEARCH
         CASE(2)
            IF (IOR== 2 .AND.  JJ/=VT%J1)   CYCLE SEARCH
            IF (IOR==-2 .AND.  JJ/=VT%J1+1) CYCLE SEARCH
            IF (II<VT%I1+1 .OR. II>VT%I2)   CYCLE SEARCH
            IF (KK<VT%K1+1 .OR. KK>VT%K2)   CYCLE SEARCH
         CASE(3)
            IF (IOR== 3 .AND.  KK/=VT%K1)   CYCLE SEARCH
            IF (IOR==-3 .AND.  KK/=VT%K1+1) CYCLE SEARCH
            IF (II<VT%I1+1 .OR. II>VT%I2)   CYCLE SEARCH
            IF (JJ<VT%J1+1 .OR. JJ>VT%J2)   CYCLE SEARCH
      END SELECT
      IF (ACTIVATE_VENT) THEN
         TW(IW) = T
      ELSE
         TW(IW) = 1000000._EB
      ENDIF
   ENDDO SEARCH

   ! Write message to .smv file

   IF (ACTIVATE_VENT)   SV_LABEL = 'OPEN_VENT'
   IF (DEACTIVATE_VENT) SV_LABEL = 'CLOSE_VENT'

   IF (N_STRINGS+2>N_STRINGS_MAX) THEN
      CALL RE_ALLOCATE_STRINGS(NM)
      STRING => MESHES(NM)%STRING
   ENDIF
   N_STRINGS = N_STRINGS + 1
   WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
   N_STRINGS = N_STRINGS + 1
   WRITE(STRING(N_STRINGS),'(I6,F10.2)') N,T

ENDDO VENT_LOOP
 
END SUBROUTINE OPEN_AND_CLOSE
 
 
SUBROUTINE CREATE_OR_REMOVE_OBST(NM,I1,I2,J1,J2,K1,K2,CR_INDEX,OBST_INDEX)

! Create or remove the obstruction whose CELLS (not nodes) are given by I1, I2, etc.

USE GEOMETRY_FUNCTIONS, ONLY : BLOCK_CELL
INTEGER :: I1,I2,J1,J2,K1,K2,I,J,K,IE,IW,IC,ICG
INTEGER, INTENT(IN) :: NM,CR_INDEX,OBST_INDEX
LOGICAL :: CREATE,REMOVE
 
CALL POINT_TO_MESH(NM)

REMOVE = .FALSE.
CREATE = .FALSE.
IF (CR_INDEX==0) REMOVE = .TRUE.
IF (CR_INDEX==1) CREATE = .TRUE.
 
! Blank or unblank cells that make up the OBSTruction
 
CALL BLOCK_CELL(NM,I1,I2,J1,J2,K1,K2,CR_INDEX,OBST_INDEX)
 
! Process the x boundaries of the OBSTruction

DO K=K1,K2
   DO J=J1,J2

      ! Nullify and/or uncover wall cells on the x=x1 side of the OBSTruction

      IF (I1 > 1) THEN
         ICG = CELL_INDEX(I1-1,J,K)
         IW  = WALL_INDEX(ICG,1)
         IF (REMOVE) BOUNDARY_TYPE(IW) = NULL_BOUNDARY
         IF (CREATE .AND. .NOT.SOLID(ICG)) CALL GET_BOUNDARY_TYPE( IW , IJKW(5,IW) , BOUNDARY_TYPE(IW) )
      ENDIF
      IC = CELL_INDEX(I1,J,K)
      IW = WALL_INDEX(IC,-1)
      IF (IW > 0 .AND. REMOVE) CALL GET_BOUNDARY_TYPE( IW , IJKW(5,IW) , BOUNDARY_TYPE(IW) )
      IF (CREATE) BOUNDARY_TYPE(IW) = NULL_BOUNDARY

      ! Nullify and/or uncover wall cells on the x=x2 side of the OBSTruction

      IF (I2 < IBAR) THEN
         ICG = CELL_INDEX(I2+1,J,K)
         IW  = WALL_INDEX(ICG,-1)
         IF (REMOVE) BOUNDARY_TYPE(IW) = NULL_BOUNDARY
         IF (CREATE .AND. .NOT.SOLID(ICG)) CALL GET_BOUNDARY_TYPE( IW , IJKW(5,IW) , BOUNDARY_TYPE(IW) )
      ENDIF
      IC = CELL_INDEX(I2,J,K)
      IW = WALL_INDEX(IC,1)
      IF (IW > 0 .AND. REMOVE) CALL GET_BOUNDARY_TYPE( IW , IJKW(5,IW) , BOUNDARY_TYPE(IW) )
      IF (CREATE) BOUNDARY_TYPE(IW) = NULL_BOUNDARY
   ENDDO
ENDDO

! Process the y boundaries of the OBSTruction

DO K=K1,K2
   DO I=I1,I2

      ! Nullify and/or uncover wall cells on the y=y1 side of the OBSTruction

      IF (J1 > 1) THEN
         ICG = CELL_INDEX(I,J1-1,K)
         IW  = WALL_INDEX(ICG,2)
         IF (REMOVE) BOUNDARY_TYPE(IW) = NULL_BOUNDARY
         IF (CREATE .AND. .NOT.SOLID(ICG)) CALL GET_BOUNDARY_TYPE( IW , IJKW(5,IW) , BOUNDARY_TYPE(IW) )
      ENDIF
      IC = CELL_INDEX(I,J1,K)
      IW = WALL_INDEX(IC,-2)
      IF (IW > 0 .AND. REMOVE) CALL GET_BOUNDARY_TYPE( IW , IJKW(5,IW) , BOUNDARY_TYPE(IW) )
      IF (CREATE) BOUNDARY_TYPE(IW) = NULL_BOUNDARY

      ! Nullify and/or uncover wall cells on the y=y2 side of the OBSTruction

      IF (J2 < JBAR) THEN
         ICG = CELL_INDEX(I,J2+1,K)
         IW  = WALL_INDEX(I