divg.f90

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

MODULE DIVG              
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
 
IMPLICIT NONE
PRIVATE
CHARACTER(255), PARAMETER :: divgid='$Id: divg.f90 645 2007-09-19 20:15:22Z drjfloyd $'
CHARACTER(255), PARAMETER :: divgrev='$Revision: 645 $'
CHARACTER(255), PARAMETER :: divgdate='$Date: 2007-09-19 16:15:22 -0400 (Wed, 19 Sep 2007) $'

PUBLIC DIVERGENCE_PART_1,DIVERGENCE_PART_2,CHECK_DIVERGENCE,GET_REV_divg
 
CONTAINS
 
 
SUBROUTINE DIVERGENCE_PART_1(T,NM)
USE COMP_FUNCTIONS, ONLY: SECOND 
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE PHYSICAL_FUNCTIONS, ONLY: GET_CP2,GET_D2,GET_K2

! Compute contributions to the divergence term
 
INTEGER, INTENT(IN) :: NM
REAL(EB), POINTER, DIMENSION(:,:,:) :: KDTDX,KDTDY,KDTDZ,DP,KP, &
          RHO_D_DYDX,RHO_D_DYDY,RHO_D_DYDZ,RHO_D,RHOP,H_RHO_D_DYDX,H_RHO_D_DYDY,H_RHO_D_DYDZ,RTRM
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP
REAL(EB) :: DELKDELT,VC,DTDX,DTDY,DTDZ,K_SUM,CP_SUM,XIF,YIF,ZIF, TNOW, &
            HDIFF,DYDX,DYDY,DYDZ,T,RDT,RHO_D_DYDN,TSI,VDOT_LEAK,TIME_RAMP_FACTOR,ZONE_VOLUME,CP_MF,Z_2,DELTA_P,PRES_RAMP_FACTOR
TYPE(SURFACE_TYPE), POINTER :: SF
INTEGER :: IW,N,IOR,II,JJ,KK,IIG,JJG,KKG,ITMP,IBC,IZ,NM_IN,II_IN,JJ_IN,KK_IN,I,J,K,IPZ,IOPZ
 
IF (SOLID_PHASE_ONLY) RETURN

TNOW=SECOND()
CALL POINT_TO_MESH(NM)
 
RDT = 1._EB/DT
 
SELECT CASE(PREDICTOR)
   CASE(.TRUE.)  
      DP => DS   
      R_PBAR = 1._EB/PBAR
      RHOP => RHOS
   CASE(.FALSE.) 
      DP => DDDT 
      R_PBAR = 1._EB/PBAR_S
      RHOP => RHO
END SELECT
 
IF (N_SPECIES > 0) THEN
   RHO_D_DYDX  => WORK1
   RHO_D_DYDY  => WORK2
   RHO_D_DYDZ  => WORK3
   SELECT CASE(PREDICTOR)
      CASE(.TRUE.)  
         YYP => YYS 
      CASE(.FALSE.) 
         YYP => YY  
   END SELECT
ENDIF

! Zero out divergence to start
 
DP  = 0._EB
IF (N_SPECIES > 0) DEL_RHO_D_DEL_Y = 0._EB

! Add species diffusion terms to divergence expression and compute diffusion term for species equations
 
SPECIES_LOOP: DO N=1,N_SPECIES
 
   ! Compute rho*D
 
   RHO_D => WORK4
    
   IF (DNS .AND. SPECIES(N)%MODE==GAS_SPECIES) THEN
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               ITMP = 0.1_EB*TMP(I,J,K)
               RHO_D(I,J,K) = RHOP(I,J,K)*SPECIES(N)%D(ITMP)
            ENDDO 
         ENDDO
      ENDDO
   ENDIF
    
   IF (DNS .AND. SPECIES(N)%MODE==MIXTURE_FRACTION_SPECIES) THEN
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               ITMP = 0.1_EB*TMP(I,J,K)
               IZ   = NINT(Z_SUM(I,J,K)*100._EB)
               IZ   = MAX(0,MIN(IZ,100))
               Z_2 = 0._EB
               IF (CO_PRODUCTION) Z_2 = YYP(I,J,K,I_PROG_CO)
               CALL GET_D2(YYP(I,J,K,I_FUEL),YYP(I,J,K,I_PROG_CO),YYP(I,J,K,I_PROG_F),Y_SUM(I,J,K),RHO_D(I,J,K),ITMP)
               RHO_D(I,J,K) = RHOP(I,J,K)*RHO_D(I,J,K)
            ENDDO
         ENDDO
      ENDDO
   ENDIF

   IF (LES) RHO_D = MU*RSC
 
   ! Compute rho*D del Y
 
   DO K=0,KBAR
      DO J=0,JBAR
         DO I=0,IBAR
            DYDX = (YYP(I+1,J,K,N)-YYP(I,J,K,N))*RDXN(I)
            RHO_D_DYDX(I,J,K) = .5_EB*(RHO_D(I+1,J,K)+RHO_D(I,J,K))*DYDX
            DYDY = (YYP(I,J+1,K,N)-YYP(I,J,K,N))*RDYN(J)
            RHO_D_DYDY(I,J,K) = .5_EB*(RHO_D(I,J+1,K)+RHO_D(I,J,K))*DYDY
            DYDZ = (YYP(I,J,K+1,N)-YYP(I,J,K,N))*RDZN(K)
            RHO_D_DYDZ(I,J,K) = .5_EB*(RHO_D(I,J,K+1)+RHO_D(I,J,K))*DYDZ
         ENDDO
      ENDDO
   ENDDO
 
   ! Correct rho*D del Y at boundaries and store rho*D at boundaries
 
   WALL_LOOP: DO IW=1,NWC
      IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_LOOP
      IIG = IJKW(6,IW) 
      JJG = IJKW(7,IW) 
      KKG = IJKW(8,IW) 
      RHODW(IW,N) = RHO_D(IIG,JJG,KKG)
      RHO_D_DYDN  = RHODW(IW,N)*(YYP(IIG,JJG,KKG,N)-YY_W(IW,N))*RDN(IW)
      DEL_RHO_D_DEL_Y(IIG,JJG,KKG,N) = DEL_RHO_D_DEL_Y(IIG,JJG,KKG,N) - RHO_D_DYDN*RDN(IW)
      IOR = IJKW(4,IW)
      SELECT CASE(IOR) 
         CASE( 1)
            RHO_D_DYDX(IIG-1,JJG,KKG) = 0._EB
         CASE(-1) 
            RHO_D_DYDX(IIG,JJG,KKG)   = 0._EB
         CASE( 2)
            RHO_D_DYDY(IIG,JJG-1,KKG) = 0._EB
         CASE(-2) 
            RHO_D_DYDY(IIG,JJG,KKG)   = 0._EB
         CASE( 3) 
            RHO_D_DYDZ(IIG,JJG,KKG-1) = 0._EB
         CASE(-3) 
            RHO_D_DYDZ(IIG,JJG,KKG)   = 0._EB
      END SELECT
   ENDDO WALL_LOOP

   ! Compute del dot h_n*rho*D del Y_n only for non-mixture fraction cases
 
   SPECIES_DIFFUSION: IF (.NOT.MIXTURE_FRACTION) THEN

      H_RHO_D_DYDX => WORK4
      H_RHO_D_DYDY => WORK5
      H_RHO_D_DYDZ => WORK6
    
      DO K=0,KBAR
         DO J=0,JBAR
            DO I=0,IBAR
               ITMP = .05_EB*(TMP(I+1,J,K)+TMP(I,J,K))
               HDIFF = SPECIES(N)%H_G(ITMP)-SPECIES(0)%H_G(ITMP)
               H_RHO_D_DYDX(I,J,K) = HDIFF*RHO_D_DYDX(I,J,K)
               ITMP = .05_EB*(TMP(I,J+1,K)+TMP(I,J,K))
               HDIFF = SPECIES(N)%H_G(ITMP)-SPECIES(0)%H_G(ITMP)
               H_RHO_D_DYDY(I,J,K) = HDIFF*RHO_D_DYDY(I,J,K)
               ITMP = .05_EB*(TMP(I,J,K+1)+TMP(I,J,K))
               HDIFF = SPECIES(N)%H_G(ITMP)-SPECIES(0)%H_G(ITMP)
               H_RHO_D_DYDZ(I,J,K) = HDIFF*RHO_D_DYDZ(I,J,K)
            ENDDO
         ENDDO
      ENDDO

      WALL_LOOP2: DO IW=1,NWC
         IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_LOOP2
         IIG = IJKW(6,IW)
         JJG = IJKW(7,IW)
         KKG = IJKW(8,IW)
         IOR  = IJKW(4,IW)
         ITMP = 0.1_EB*TMP_F(IW)
         HDIFF = SPECIES(N)%H_G(ITMP)-SPECIES(0)%H_G(ITMP)
         RHO_D_DYDN = RHODW(IW,N)*(YYP(IIG,JJG,KKG,N)-YY_W(IW,N))*RDN(IW)
         SELECT CASE(IOR)
            CASE( 1) 
               H_RHO_D_DYDX(IIG-1,JJG,KKG) = 0._EB
            CASE(-1) 
               H_RHO_D_DYDX(IIG,JJG,KKG)   = 0._EB
            CASE( 2) 
               H_RHO_D_DYDY(IIG,JJG-1,KKG) = 0._EB
            CASE(-2) 
               H_RHO_D_DYDY(IIG,JJG,KKG)   = 0._EB
            CASE( 3) 
               H_RHO_D_DYDZ(IIG,JJG,KKG-1) = 0._EB
            CASE(-3) 
               H_RHO_D_DYDZ(IIG,JJG,KKG)   = 0._EB
         END SELECT
         DP(IIG,JJG,KKG) = DP(IIG,JJG,KKG) - RDN(IW)*HDIFF*RHO_D_DYDN
      ENDDO WALL_LOOP2
 
      CYLINDER: SELECT CASE(CYLINDRICAL)
         CASE(.FALSE.) CYLINDER  ! 3D or 2D Cartesian Coords
            DO K=1,KBAR
               DO J=1,JBAR
                  DO I=1,IBAR
                     DP(I,J,K) = DP(I,J,K) + & 
                                 (H_RHO_D_DYDX(I,J,K)-H_RHO_D_DYDX(I-1,J,K))*RDX(I) + &
                                 (H_RHO_D_DYDY(I,J,K)-H_RHO_D_DYDY(I,J-1,K))*RDY(J) + &
                                 (H_RHO_D_DYDZ(I,J,K)-H_RHO_D_DYDZ(I,J,K-1))*RDZ(K)
                  ENDDO
               ENDDO
            ENDDO
         CASE(.TRUE.) CYLINDER  ! 2D Cylindrical Coords
            J = 1
            DO K=1,KBAR
               DO I=1,IBAR
                  DP(I,J,K) = DP(I,J,K) + &
                  (R(I)*H_RHO_D_DYDX(I,J,K)-R(I-1)*H_RHO_D_DYDX(I-1,J,K))*RDX(I)*RRN(I) + &
                       (H_RHO_D_DYDZ(I,J,K)-       H_RHO_D_DYDZ(I,J,K-1))*RDZ(K)
               ENDDO
            ENDDO
      END SELECT CYLINDER
    
   ENDIF SPECIES_DIFFUSION

   ! Compute del dot rho*D del Y_n or del dot rho D del Z
 
   CYLINDER2: SELECT CASE(CYLINDRICAL)
      CASE(.FALSE.) CYLINDER2  ! 3D or 2D Cartesian Coords
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  DEL_RHO_D_DEL_Y(I,J,K,N) = DEL_RHO_D_DEL_Y(I,J,K,N) + &
                                 (RHO_D_DYDX(I,J,K)-RHO_D_DYDX(I-1,J,K))*RDX(I) + &
                                 (RHO_D_DYDY(I,J,K)-RHO_D_DYDY(I,J-1,K))*RDY(J) + &
                                 (RHO_D_DYDZ(I,J,K)-RHO_D_DYDZ(I,J,K-1))*RDZ(K)
               ENDDO
            ENDDO 
         ENDDO
      CASE(.TRUE.) CYLINDER2  ! 2D Cylindrical Coords
         J=1
         DO K=1,KBAR
            DO I=1,IBAR
               DEL_RHO_D_DEL_Y(I,J,K,N) = DEL_RHO_D_DEL_Y(I,J,K,N) + &
               (R(I)*RHO_D_DYDX(I,J,K)-R(I-1)*RHO_D_DYDX(I-1,J,K))*RDX(I)*RRN(I) + &
                    (RHO_D_DYDZ(I,J,K)-       RHO_D_DYDZ(I,J,K-1))*RDZ(K)
            ENDDO
         ENDDO
   END SELECT CYLINDER2

   ! Compute -Sum h_n del dot rho*D del Y_n
 
   SPECIES_DIFFUSION_2: IF (.NOT.MIXTURE_FRACTION) THEN
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               ITMP = 0.1_EB*TMP(I,J,K)
               HDIFF = SPECIES(N)%H_G(ITMP)-SPECIES(0)%H_G(ITMP)
               DP(I,J,K) = DP(I,J,K) - HDIFF*DEL_RHO_D_DEL_Y(I,J,K,N)
            ENDDO
         ENDDO
      ENDDO
   ENDIF SPECIES_DIFFUSION_2
 
ENDDO SPECIES_LOOP
 

! Compute del dot k del T
 
ENERGY: IF (.NOT.ISOTHERMAL) THEN
 
   KDTDX => WORK1
   KDTDY => WORK2
   KDTDZ => WORK3
   KP    => WORK4
 
   ! Compute thermal conductivity k (KP)
 
   K_DNS_OR_LES: IF (DNS) THEN
 
      IF (.NOT.MIXTURE_FRACTION) THEN
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  ITMP = 0.1_EB*TMP(I,J,K)
                  K_SUM = SPECIES(0)%K(ITMP)
                  DO N=1,N_SPECIES
                     K_SUM = K_SUM + YYP(I,J,K,N)*(SPECIES(N)%K(ITMP)-SPECIES(0)%K(ITMP))
                  END DO
                  KP(I,J,K) = K_SUM
               ENDDO
            ENDDO
         ENDDO
      ELSE   ! MIXTURE_FRACTION TRUE
         Z_2 = 0._EB
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  ITMP = 0.1_EB*TMP(I,J,K)
                  IF(CO_PRODUCTION) THEN
                     CALL GET_CP2(YY(I,J,K,I_FUEL),YY(I,J,K,I_PROG_CO),YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),CP_MF,ITMP)
                     IF (N_SPECIES > 3) THEN
                        CP_SUM = 0._EB
                        DO N=1,N_SPECIES
                           IF (SPECIES(N)%MODE/=MIXTURE_FRACTION_SPECIES) &
                           CP_SUM = CP_SUM + YYP(I,J,K,N)*(SPECIES(N)%CP(ITMP)-SPECIES(0)%CP(ITMP))
                        END DO
                        CP_MF = CP_SUM+(1._EB-Y_SUM(I,J,K))*CP_MF
                     ENDIF
                  ELSE
                     CALL GET_CP2(YY(I,J,K,I_FUEL),0._EB,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),CP_MF,ITMP)                  
                     IF (N_SPECIES > 2) THEN
                        CP_SUM = 0._EB
                        DO N=1,N_SPECIES
                           IF (SPECIES(N)%MODE/=MIXTURE_FRACTION_SPECIES) &
                           CP_SUM = CP_SUM + YYP(I,J,K,N)*(SPECIES(N)%CP(ITMP)-SPECIES(0)%CP(ITMP))
                        END DO
                        CP_MF = CP_SUM+(1._EB-Y_SUM(I,J,K))*CP_MF
                     ENDIF
                  ENDIF
                  KP(I,J,K) = MU(I,J,K)*CP_MF*RPR
               ENDDO
            ENDDO
         ENDDO
      ENDIF
      BOUNDARY_LOOP: DO IW=1,NEWC
         II  = IJKW(1,IW)
         JJ  = IJKW(2,IW)
         KK  = IJKW(3,IW)
         IIG = IJKW(6,IW)
         JJG = IJKW(7,IW)
         KKG = IJKW(8,IW)
         KP(II,JJ,KK) = KP(IIG,JJG,KKG)
      ENDDO BOUNDARY_LOOP
 
   ELSE K_DNS_OR_LES
    
      KP = MU*CPOPR
 
   ENDIF K_DNS_OR_LES
 
   ! Compute k*dT/dx, etc

   DO K=0,KBAR
      DO J=0,JBAR
         DO I=0,IBAR
            DTDX = (TMP(I+1,J,K)-TMP(I,J,K))*RDXN(I)
            KDTDX(I,J,K) = .5_EB*(KP(I+1,J,K)+KP(I,J,K))*DTDX
            DTDY = (TMP(I,J+1,K)-TMP(I,J,K))*RDYN(J)
            KDTDY(I,J,K) = .5_EB*(KP(I,J+1,K)+KP(I,J,K))*DTDY
            DTDZ = (TMP(I,J,K+1)-TMP(I,J,K))*RDZN(K)
            KDTDZ(I,J,K) = .5_EB*(KP(I,J,K+1)+KP(I,J,K))*DTDZ
         ENDDO
      ENDDO
   ENDDO

   ! Correct thermal gradient (k dT/dn) at boundaries
 
   CORRECTION_LOOP: DO IW=1,NWC
      IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE CORRECTION_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) 
      KW(IW) = KP(IIG,JJG,KKG)
      IOR = IJKW(4,IW)
      SELECT CASE(IOR)
         CASE( 1)
            KDTDX(II,JJ,KK)   = 0._EB
         CASE(-1)
            KDTDX(II-1,JJ,KK) = 0._EB
         CASE( 2)
            KDTDY(II,JJ,KK)   = 0._EB
         CASE(-2)
            KDTDY(II,JJ-1,KK) = 0._EB
         CASE( 3)
            KDTDZ(II,JJ,KK)   = 0._EB
         CASE(-3)
            KDTDZ(II,JJ,KK-1) = 0._EB
      END SELECT
      DP(IIG,JJG,KKG) = DP(IIG,JJG,KKG) - QCONF(IW)*RDN(IW)
   ENDDO CORRECTION_LOOP

   ! Compute (q + del dot k del T) and add to the divergence
 
   CYLINDER3: SELECT CASE(CYLINDRICAL)
      CASE(.FALSE.) CYLINDER3   ! 3D or 2D Cartesian
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  DELKDELT = (KDTDX(I,J,K)-KDTDX(I-1,J,K))*RDX(I) + &
                             (KDTDY(I,J,K)-KDTDY(I,J-1,K))*RDY(J) + &
                             (KDTDZ(I,J,K)-KDTDZ(I,J,K-1))*RDZ(K)
                  DP(I,J,K) = DP(I,J,K) + DELKDELT + Q(I,J,K) + QR(I,J,K)
               ENDDO 
            ENDDO
         ENDDO
      CASE(.TRUE.) CYLINDER3   ! 2D Cylindrical
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  DELKDELT = & 
                  (R(I)*KDTDX(I,J,K)-R(I-1)*KDTDX(I-1,J,K))*RDX(I)*RRN(I) + &
                       (KDTDZ(I,J,K)-       KDTDZ(I,J,K-1))*RDZ(K)
                  DP(I,J,K) = DP(I,J,K) + DELKDELT + Q(I,J,K) + QR(I,J,K)
               ENDDO 
            ENDDO
         ENDDO
   END SELECT CYLINDER3
 
ENDIF ENERGY

! Compute RTRM = R*sum(Y_i/M_i)/(PBAR*C_P) and multiply it by divergence terms already summed up
 
RTRM => WORK1
 
IF (MIXTURE_FRACTION) THEN
   Z_2 = 0._EB
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            ITMP = 0.1_EB*TMP(I,J,K)
            IF(CO_PRODUCTION) THEN
               CALL GET_CP2(YY(I,J,K,I_FUEL),YY(I,J,K,I_PROG_CO),YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),CP_MF,ITMP)
               IF (N_SPECIES > 3) THEN
                  CP_SUM = 0._EB
                  DO N=1,N_SPECIES
                     IF (SPECIES(N)%MODE/=MIXTURE_FRACTION_SPECIES) &
                     CP_SUM = CP_SUM + YYP(I,J,K,N)*(SPECIES(N)%CP(ITMP)-SPECIES(0)%CP(ITMP))
                     CP_MF = CP_SUM+(1._EB-Y_SUM(I,J,K))*CP_MF
                  END DO
               ENDIF
            ELSE
               CALL GET_CP2(YY(I,J,K,I_FUEL),Z_2,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),CP_MF,ITMP)  
               IF (N_SPECIES > 2) THEN