wall.f90

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

MODULE WALL_ROUTINES
 
! Compute the wall boundary conditions
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
 
IMPLICIT NONE
PRIVATE
CHARACTER(255), PARAMETER :: wallid='$Id: wall.f90 721 2007-10-01 18:33:34Z drjfloyd $'
CHARACTER(255), PARAMETER :: wallrev='$Revision: 721 $'
CHARACTER(255), PARAMETER :: walldate='$Date: 2007-10-01 14:33:34 -0400 (Mon, 01 Oct 2007) $'

PUBLIC WALL_BC,GET_REV_wall
 
CONTAINS

SUBROUTINE WALL_BC(T,NM)
USE COMP_FUNCTIONS, ONLY: SECOND
REAL(EB) :: TNOW
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM

IF (EVACUATION_ONLY(NM)) RETURN

TNOW=SECOND()

CALL POINT_TO_MESH(NM)

IF (.NOT. ISOTHERMAL) CALL THERMAL_BC(T)
CALL SPECIES_BC(T)
CALL DENSITY_BC

TUSED(6,NM)=TUSED(6,NM)+SECOND()-TNOW


CONTAINS
 

SUBROUTINE THERMAL_BC(T)

! Thermal boundary conditions for adiabatic, fixed temperature, fixed flux and interpolated boundaries.
! One dimensional heat transfer and pyrolysis is done in PYROLYSIS, which is called at the end of this routine.

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP 
REAL(EB) :: DT_BC,T,TSI,TMP_G,DTMP,TMP_OTHER,CP_TERM,RHOWAL,RAMP_FACTOR,QNET,FDERIV,TMP_EXTERIOR
INTEGER  :: IOR,II,JJ,KK,IBC,IIG,JJG,KKG, IW
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (VENTS_TYPE), POINTER :: VT
 
IF (PREDICTOR) THEN
   UU => U
   VV => V
   WW => W
   RHOP => RHOS
ELSE
   UU => US
   VV => VS
   WW => WS
   RHOP => RHO
ENDIF
 
! Loop through all boundary cells and apply heat transfer method, except for thermally-thick cells
 
HEAT_FLUX_LOOP: DO IW=1,NWC
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE HEAT_FLUX_LOOP
   II  = IJKW(1,IW)
   JJ  = IJKW(2,IW)
   KK  = IJKW(3,IW)
   IIG = IJKW(6,IW)
   JJG = IJKW(7,IW)
   KKG = IJKW(8,IW)
   IOR = IJKW(4,IW)
   IBC = IJKW(5,IW)
   SF  => SURFACE(IBC)
 
   METHOD_OF_HEAT_TRANSFER: SELECT CASE(SF%THERMAL_BC_INDEX)
 
      CASE (NO_CONVECTION) METHOD_OF_HEAT_TRANSFER

         TMP_F(IW) = TMP(IIG,JJG,KKG)
         TMP_W(IW) = TMP_F(IW)

      CASE (INFLOW_OUTFLOW) METHOD_OF_HEAT_TRANSFER 
 
         TMP_F(IW) = TMP(IIG,JJG,KKG)
         TMP_EXTERIOR = TMP_0(KK)
         IF (VENT_INDEX(IW)>0) THEN
            VT => VENTS(VENT_INDEX(IW))
            IF (VT%TMP_EXTERIOR>0._EB) TMP_EXTERIOR = VT%TMP_EXTERIOR
         ENDIF
         SELECT CASE(IOR)
            CASE( 1) 
               IF (UU(II,JJ,KK)>=0._EB)   TMP_F(IW) = TMP_EXTERIOR
            CASE(-1) 
               IF (UU(II-1,JJ,KK)<=0._EB) TMP_F(IW) = TMP_EXTERIOR
            CASE( 2) 
               IF (VV(II,JJ,KK)>=0._EB)   TMP_F(IW) = TMP_EXTERIOR
            CASE(-2) 
               IF (VV(II,JJ-1,KK)<=0._EB) TMP_F(IW) = TMP_EXTERIOR
            CASE( 3) 
               IF (WW(II,JJ,KK)>=0._EB)   TMP_F(IW) = TMP_EXTERIOR
            CASE(-3) 
               IF (WW(II,JJ,KK-1)<=0._EB) TMP_F(IW) = TMP_EXTERIOR
            END SELECT
         TMP(II,JJ,KK) = TMP_F(IW)
         TMP_W(IW) = TMP_F(IW)
 
      CASE (SPECIFIED_TEMPERATURE) METHOD_OF_HEAT_TRANSFER

         TMP_G = TMP(IIG,JJG,KKG)
         IF (TW(IW)==T_BEGIN) THEN
            TSI = T
         ELSE
            TSI = T - TW(IW)
         ENDIF
         TMP_F(IW) = TMP_0(KK) + EVALUATE_RAMP(TSI,SF%TAU(TIME_TEMP),SF%RAMP_INDEX(TIME_TEMP))*(SF%TMP_FRONT-TMP_0(KK))
         DTMP = TMP_G - TMP_F(IW)
         HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)
         QCONF(IW) = HEAT_TRANS_COEF(IW)*DTMP
         RHOWAL    = 0.5_EB*(RHOP(IIG,JJG,KKG)+RHO_W(IW))
         CP_TERM   = MAX(0._EB,-CP_GAMMA*UW(IW)*RHOWAL)
         TMP_W(IW) = ( (RDN(IW)*KW(IW)-0.5_EB*CP_TERM)*TMP_G + CP_TERM*TMP_F(IW)-QCONF(IW) )/(0.5_EB*CP_TERM+RDN(IW)*KW(IW))
         TMP_W(IW) = MAX(TMPMIN,TMP_W(IW))

      CASE (ADIABATIC_INDEX) METHOD_OF_HEAT_TRANSFER

         TMP_G = TMP(IIG,JJG,KKG)
         TMP_OTHER = TMP_F(IW)
         ADLOOP: DO
            DTMP = TMP_G - TMP_OTHER
            HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)
            IF (RADIATION) THEN
               QNET = HEAT_TRANS_COEF(IW)*DTMP + QRADIN(IW) - E_WALL(IW) * SIGMA * TMP_OTHER ** 4
               FDERIV = -HEAT_TRANS_COEF(IW) -  4._EB * E_WALL(IW) * SIGMA * TMP_OTHER ** 3
            ELSE
               QNET = HEAT_TRANS_COEF(IW)*DTMP
               FDERIV = -HEAT_TRANS_COEF(IW)
            ENDIF
            IF (FDERIV /= 0._EB) TMP_OTHER = TMP_OTHER - QNET / FDERIV
            IF (ABS(TMP_OTHER - TMP_F(IW)) / TMP_F(IW) < 0.0001) THEN
               TMP_F(IW) = TMP_OTHER
               EXIT ADLOOP
            ELSE
               TMP_F(IW) = TMP_OTHER
               CYCLE ADLOOP
            ENDIF
         ENDDO ADLOOP
         DTMP = TMP_G - TMP_F(IW)
         HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)
         QCONF(IW) = HEAT_TRANS_COEF(IW)*DTMP
         RHOWAL    = 0.5_EB*(RHOP(IIG,JJG,KKG)+RHO_W(IW))
         CP_TERM   = MAX(0._EB,-CP_GAMMA*UW(IW)*RHOWAL)
         TMP_W(IW) = ( (RDN(IW)*KW(IW)-0.5_EB*CP_TERM)*TMP_G + CP_TERM*TMP_F(IW)-QCONF(IW) )/(0.5_EB*CP_TERM+RDN(IW)*KW(IW))
         TMP_W(IW) = MAX(TMPMIN,TMP_W(IW))

      CASE (SPECIFIED_HEAT_FLUX) METHOD_OF_HEAT_TRANSFER
         IF (TW(IW)==T_BEGIN) THEN
            TSI = T
         ELSE
            TSI = T - TW(IW)
         ENDIF
         RAMP_FACTOR = EVALUATE_RAMP(TSI,SF%TAU(TIME_HEAT),SF%RAMP_INDEX(TIME_HEAT))
         TMP_F(IW) = TMPA + RAMP_FACTOR*(SF%TMP_FRONT-TMPA)
         QCONF(IW) = -RAMP_FACTOR*SF%CONVECTIVE_HEAT_FLUX*AREA_ADJUST(IW)
         RHOWAL    = 0.5_EB*(RHOP(IIG,JJG,KKG)+RHO_W(IW))
         TMP_G     = TMP(IIG,JJG,KKG)
         CP_TERM   = MAX(0._EB,-CP_GAMMA*UW(IW)*RHOWAL)
         TMP_W(IW) = ( (RDN(IW)*KW(IW)-0.5_EB*CP_TERM)*TMP_G + CP_TERM*TMP_F(IW)-QCONF(IW) )/(0.5_EB*CP_TERM+RDN(IW)*KW(IW))
         TMP_W(IW) = MAX(TMPMIN,TMP_W(IW))
 
      CASE (INTERPOLATED_BC) METHOD_OF_HEAT_TRANSFER
 
         TMP_G = TMP(IIG,JJG,KKG)
         TMP_OTHER =OMESH(IJKW(9,IW))%TMP(IJKW(10,IW),IJKW(11,IW),IJKW(12,IW))
         IF (CELL_VOLUME_RATIO(IW)<0.5 .OR.  CELL_VOLUME_RATIO(IW)>2.0) THEN
            TMP_W(IW) = TMP_G
            SELECT CASE(IOR)
               CASE( 1) 
                  IF (UU(II,JJ,KK)>=0._EB)   TMP_W(IW) = TMP_OTHER
               CASE(-1) 
                  IF (UU(II-1,JJ,KK)<=0._EB) TMP_W(IW) = TMP_OTHER
               CASE( 2) 
                  IF (VV(II,JJ,KK)>=0._EB)   TMP_W(IW) = TMP_OTHER
               CASE(-2) 
                  IF (VV(II,JJ-1,KK)<=0._EB) TMP_W(IW) = TMP_OTHER
               CASE( 3) 
                  IF (WW(II,JJ,KK)>=0._EB)   TMP_W(IW) = TMP_OTHER
               CASE(-3) 
                  IF (WW(II,JJ,KK-1)<=0._EB) TMP_W(IW) = TMP_OTHER
            END SELECT
         ELSE
            TMP_W(IW) = TMP_OTHER
         ENDIF
         TMP_F(IW) = TMP_W(IW)
         QCONF(IW) = KW(IW)*(TMP_G-TMP_W(IW))*RDN(IW)
         TMP(II,JJ,KK) = TMP_W(IW)
 
   END SELECT METHOD_OF_HEAT_TRANSFER
 
! Record wall temperature in the ghost cell
 
   IF (SOLID(CELL_INDEX(II,JJ,KK))) TMP(II,JJ,KK) = MAX(100._EB,MIN(4900._EB,TMP_W(IW)))
 
ENDDO HEAT_FLUX_LOOP
 
! For thermally-thick boundary conditions, call the routine PYROLYSIS
 
IF (CORRECTOR) THEN
   WALL_COUNTER = WALL_COUNTER + 1
   IF (WALL_COUNTER==WALL_INCREMENT) THEN
      DT_BC    = T - BC_CLOCK
      BC_CLOCK = T
      CALL PYROLYSIS(T,DT_BC)
      WALL_COUNTER = 0
   ENDIF
ENDIF
 
END SUBROUTINE THERMAL_BC
 
 
SUBROUTINE SPECIES_BC(T)
USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP
 
REAL(EB) T,YY_WALL,YY_G,DENOM,YY_OTHER(MAX_SPECIES),YY_EXTERIOR(MAX_SPECIES),RHO_G,UN,DD,EPSB,MFT,TSI
INTEGER IBC,IIG,JJG,KKG,IOR,IWB,IW, II, JJ, KK, N
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (VENTS_TYPE), POINTER :: VT
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP
 
IF (PREDICTOR) THEN
   UU => U
   VV => V
   WW => W
   RHOP => RHOS
   IF (N_SPECIES > 0) YYP => YYS
ELSE
   UU => US
   VV => VS
   WW => WS
   RHOP => RHO
   IF (N_SPECIES > 0) YYP => YY
ENDIF 
 
! Loop through the wall cells, apply mass boundary conditions

WALL_CELL_LOOP: DO IW=1,NWC
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_CELL_LOOP
   IBC = IJKW(5,IW)
   SF  => SURFACE(IBC)
   IF (N_SPECIES==0 .AND. .NOT. SF%SPECIES_BC_INDEX==SPECIFIED_MASS_FLUX) CYCLE WALL_CELL_LOOP
   IF (N_SPECIES==0 .AND. SF%SPECIES_BC_INDEX==SPECIFIED_MASS_FLUX .AND. SF%MASS_FLUX(0) == 0._EB) CYCLE WALL_CELL_LOOP
   II  = IJKW(1,IW)
   JJ  = IJKW(2,IW)
   KK  = IJKW(3,IW)
   IOR = IJKW(4,IW)
   IIG = IJKW(6,IW)
   JJG = IJKW(7,IW)
   KKG = IJKW(8,IW)
   
   METHOD_OF_MASS_TRANSFER: SELECT CASE(SF%SPECIES_BC_INDEX)
 
      CASE (NO_MASS_FLUX) METHOD_OF_MASS_TRANSFER 

         IF (.NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) YY_W(IW,1:N_SPECIES) = YYP(IIG,JJG,KKG,1:N_SPECIES)
         IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
            YY_EXTERIOR(1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
            VT => VENTS(VENT_INDEX(IW))
            DO N=1,N_SPECIES
               IF (VT%MASS_FRACTION(N)>-1._EB) YY_EXTERIOR(N) = VT%MASS_FRACTION(N)
            ENDDO
            SELECT CASE(IOR)
               CASE( 1) 
                  IF (UU(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)
               CASE(-1)
                  IF (UU(II-1,JJ,KK)<0._EB) YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)
               CASE( 2) 
                  IF (VV(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)
               CASE(-2)
                  IF (VV(II,JJ-1,KK)<0._EB) YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)
               CASE( 3)
                  IF (WW(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)
               CASE(-3)
                  IF (WW(II,JJ,KK-1)<0._EB) YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)
            END SELECT
            YYP(II,JJ,KK,1:N_SPECIES)=YY_W(IW,1:N_SPECIES)
         ENDIF
         IF ( SF%LEAK_PATH(1)>-1 .AND. UWS(IW)<0._EB) THEN
            RHO_G = RHOP(IIG,JJG,KKG)
            UN = -UWS(IW)
            IF (PREDICTOR) EPSB = -.5_EB*UN**2*DT*RDN(IW)
            IF (CORRECTOR) EPSB =  .5_EB*UN**2*DT*RDN(IW)
            SPECIES_LOOP_1: DO N=1,N_SPECIES
               DD    = RHODW(IW,N)*RDN(IW)
               YY_G  = YYP(IIG,JJG,KKG,N)
               DENOM = DD + (.5_EB*UN+EPSB)*RHO_W(IW)
               YY_W(IW,N) = ( MASSFLUX(IW,N) + YY_G*(DD + (EPSB-.5_EB*UN)*RHO_G) ) / DENOM
            ENDDO SPECIES_LOOP_1
         ENDIF
 
      CASE (SPECIFIED_MASS_FRACTION) METHOD_OF_MASS_TRANSFER
         IF (TW(IW)==T_BEGIN) THEN
            IF (PREDICTOR) TSI = T + DT
            IF (CORRECTOR) TSI = T
         ELSE
            IF (PREDICTOR) TSI = T + DT - TW(IW)
            IF (CORRECTOR) TSI = T      - TW(IW)
         ENDIF
         DO N=1,N_SPECIES
            IF (UWS(IW) <= 0._EB) THEN
               YY_WALL = SPECIES(N)%YY0 + EVALUATE_RAMP(TSI,SF%TAU(N),SF%RAMP_INDEX(N))*(SF%MASS_FRACTION(N)-SPECIES(N)%YY0)
            ELSE
               YY_WALL = YYP(IIG,JJG,KKG,N)
            ENDIF
            IF (DNS) YY_W(IW,N) = 2._EB*YY_WALL - YYP(IIG,JJG,KKG,N)
            IF (LES) YY_W(IW,N) =       YY_WALL 
         ENDDO
 
      CASE (SPECIFIED_MASS_FLUX) METHOD_OF_MASS_TRANSFER

         ! If the current time is before the "activation" time, TW, apply simple BCs and get out

         IF (T < TW(IW) .AND. N_SPECIES > 0) THEN
            IF (.NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) YY_W(IW,1:N_SPECIES)      = YYP(IIG,JJG,KKG,1:N_SPECIES)
            IF (SOLID(CELL_INDEX(II,JJ,KK)))         YYP(II,JJ,KK,1:N_SPECIES) = YY_W(IW,1:N_SPECIES)
            IF (PREDICTOR) UWS(IW)   = 0._EB 
            MASSFLUX(IW,1:N_SPECIES) = 0._EB
            ACTUAL_BURN_RATE(IW)     = 0._EB
            CYCLE WALL_CELL_LOOP
         ENDIF

         MFT = 0._EB

         ! If the user has specified the burning rate, evaluate the ramp and other related parameters

         SUM_MASSFLUX_LOOP: DO N=0,N_SPECIES
            IF (SF%MASS_FLUX(N) > 0._EB) THEN  ! Use user-specified ramp-up of mass flux
               IF (TW(IW)==T_BEGIN) THEN
                  TSI = T
               ELSE
                  TSI = T - TW(IW)
               ENDIF
               MASSFLUX(IW,N) = EVALUATE_RAMP(TSI,SF%TAU(N),SF%RAMP_INDEX(N))*SF%MASS_FLUX(N)
               IF (N==I_FUEL) THEN
                  IF (EW(IW)>0._EB) MASSFLUX(IW,N) = MASSFLUX(IW,N)*EXP(-EW(IW))
                  ACTUAL_BURN_RATE(IW) = MASSFLUX(IW,N)
               ENDIF
            ENDIF
            MASSFLUX(IW,N) = MASSFLUX(IW,N)*AREA_ADJUST(IW)
            IF (N==I_FUEL) ACTUAL_BURN_RATE(IW) = ACTUAL_BURN_RATE(IW)*AREA_ADJUST(IW)
            MFT = MFT + MASSFLUX(IW,N)
         ENDDO SUM_MASSFLUX_LOOP
 
         ! Add total fuel consumed to various summing arrays

         CONSUME_FUEL: IF (CORRECTOR .AND. SF%THERMALLY_THICK .AND. I_FUEL /= 0) THEN  
            IF (SF%SURFACE_DENSITY>0._EB .AND. SF%BACKING==EXPOSED) THEN
               IWB = WALL_INDEX_BACK(IW)
               IF (BOUNDARY_TYPE(IWB)==SOLID_BOUNDARY) MASS_LOSS(IWB) = MASS_LOSS(IWB) + MASSFLUX(IW,I_FUEL)*DT
            ENDIF
            MASS_LOSS(IW) = MASS_LOSS(IW) + MASSFLUX(IW,I_FUEL)*DT
            OBSTRUCTION(OBST_INDEX_W(IW))%MASS = OBSTRUCTION(OBST_INDEX_W(IW))%MASS - MASSFLUX(IW,I_FUEL)*DT*AW(IW)
         ENDIF CONSUME_FUEL

         ! Compute the ghost cell value of the species to get the right mass flux
 
         RHO_G = RHOP(IIG,JJG,KKG)
         UN    = 2._EB*MFT/(RHO_W(IW)+RHO_G)
         IF (PREDICTOR) UWS(IW) = -UN
         IF (PREDICTOR) EPSB = -.5_EB*UN**2*DT*RDN(IW)
         IF (CORRECTOR) EPSB =  .5_EB*UN**2*DT*RDN(IW)
         SPECIES_LOOP: DO N=1,N_SPECIES
            DD    = RHODW(IW,N)*RDN(IW)
            YY_G  = YYP(IIG,JJG,KKG,N)
            DENOM = DD + (.5_EB*UN+EPSB)*RHO_W(IW)
            YY_W(IW,N) = ( MASSFLUX(IW,N) + YY_G*(DD + (EPSB-.5_EB*UN)*RHO_G) ) / DENOM
         ENDDO SPECIES_LOOP

      CASE (INTERPOLATED_BC) METHOD_OF_MASS_TRANSFER
 
         YY_OTHER(1:N_SPECIES) = OMESH(IJKW(9,IW))%YY(IJKW(10,IW),IJKW(11,IW),IJKW(12,IW),1:N_SPECIES)
         IF (CELL_VOLUME_RATIO(IW)<0.5_EB .OR. CELL_VOLUME_RATIO(IW)>2._EB) THEN
            IF (.NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) YY_W(IW,1:N_SPECIES) = YYP(IIG,JJG,KKG,1:N_SPECIES) 
            SELECT CASE(IOR)
               CASE( 1)