velo.f90

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

ENDDO
 
! Adjust FVX, FVY and FVZ at solid, internal obstructions for no flux
 
CALL NO_FLUX
 
END SUBROUTINE VELOCITY_FLUX_ISOTHERMAL
 
 
SUBROUTINE VELOCITY_FLUX_CYLINDRICAL(T,NM)
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP 
USE PHYSICAL_FUNCTIONS, ONLY: GET_MU2
! Compute convective and diffusive terms for 2D axisymmetric
 
REAL(EB) :: T,DMUDX
INTEGER :: I0
INTEGER, INTENT(IN) :: NM
REAL(EB) :: MUY,UP,UM,WP,WM,VTRM, &
            DTXZDZ,DTXZDX, &
            DUDX,DWDZ,DUDZ,DWDX,WOMY,UOMY,PMDT,MPDT, &
            S2,C,CDXDYDZTT,AH,RRHO,GX,GZ,&
            TXXP,TXXM,TZZP,TZZM,DTXXDX,DTZZDZ, &
            EPSUP,EPSUM,EPSWP,EPSWM,MU_SUM,DUMMY=0._EB
INTEGER :: II,JJ,KK,I,J,K,IW,IIG,JJG,KKG,ITMP,N,IE,IEC
REAL(EB), POINTER, DIMENSION(:,:,:) :: TXZ,OMY,UU,WW,RHOP,DP
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP
 
CALL POINT_TO_MESH(NM)
 
IF (PREDICTOR) THEN
   UU => U
   WW => W
   DP => D  
   RHOP => RHO
   IF (N_SPECIES>0) YYP => YY
ELSE
   UU => US
   WW => WS
   DP => DS
   RHOP => RHOS
   IF (N_SPECIES>0) YYP => YYS
ENDIF
 
TXZ => WORK2
OMY => WORK5
 
CALC_MU: IF (PREDICTOR) THEN    
 
   LES_VS_DNS: IF (LES) THEN   ! Smagorinsky model (LES)

      C = CSMAG**2
      J = 1
      KLOOP: DO K=1,KBAR
         ILOOP: DO I=1,IBAR
            IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ILOOP
            CDXDYDZTT = C*DX(I)*DZ(K)
            DUDX = RDX(I)*(UU(I,J,K)-UU(I-1,J,K))
            DWDZ = RDZ(K)*(WW(I,J,K)-WW(I,J,K-1))
            DWDX = 0.25_EB*RDX(I)*(WW(I+1,J,K)-WW(I-1,J,K)+WW(I+1,J,K-1)-WW(I-1,J,K-1))
            DUDZ = 0.25_EB*RDZ(K)*(UU(I,J,K+1)-UU(I,J,K-1)+UU(I-1,J,K+1)-UU(I-1,J,K-1))
            S2   = 2._EB*(DUDX*DUDX     +  DWDZ*DWDZ ) +  (DUDZ+DWDX)**2 -  TWTH*DP(I,J,K)**2
            S2   = MAX(0._EB,S2)
            ITMP = 0.1_EB*TMP(I,J,K)
            MU(I,J,K) = MAX(SPECIES(0)%MU(ITMP), RHOP(I,J,K)*CDXDYDZTT*SQRT(S2))
         ENDDO ILOOP
      ENDDO KLOOP

   ELSE LES_VS_DNS           ! DNS viscosity

      MIXTURE_FRACTION_DNS: IF (.NOT.MIXTURE_FRACTION) THEN
         J = 1
         DO K=1,KBAR
            IVLOOP: DO I=1,IBAR
               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE IVLOOP
               ITMP = 0.1_EB*TMP(I,J,K)
               MU_SUM = SPECIES(0)%MU(ITMP)
               DO N=1,N_SPECIES
                  MU_SUM = MU_SUM + YYP(I,J,K,N)*(SPECIES(N)%MU(ITMP)-SPECIES(0)%MU(ITMP))
               ENDDO
               MU(I,J,K) = MU_SUM
            ENDDO IVLOOP
         ENDDO
      ELSE MIXTURE_FRACTION_DNS
         J = 1
         DO K=1,KBAR
            IVLOOP2: DO I=1,IBAR
               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE IVLOOP2
               ITMP = 0.1_EB*TMP(I,J,K)
               IF(CO_PRODUCTION) THEN
                  CALL GET_MU2(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),MU(I,J,K),ITMP)
               ELSE
                  CALL GET_MU2(YY(I,J,K,I_FUEL),0._EB,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),MU(I,J,K),ITMP)                  
               ENDIF
            ENDDO IVLOOP2
         ENDDO
      ENDIF MIXTURE_FRACTION_DNS
 
   ENDIF LES_VS_DNS
 
! Mirror viscosity into solids
 
   DO IW=1,NWC
      II  = IJKW(1,IW)
      JJ  = IJKW(2,IW)
      KK  = IJKW(3,IW)
      IIG = IJKW(6,IW)
      JJG = IJKW(7,IW)
      KKG = IJKW(8,IW)
      MU(II,JJ,KK) = MU(IIG,JJG,KKG)
   ENDDO
 
ENDIF CALC_MU
 
! Compute vorticity and stress tensor components
 
DO K=0,KBAR
   DO J=0,JBAR
      DO I=0,IBAR
         DUDZ = RDZN(K)*(UU(I,J,K+1)-UU(I,J,K))
         DWDX = RDXN(I)*(WW(I+1,J,K)-WW(I,J,K))
         OMY(I,J,K) = DUDZ - DWDX
         MUY = 0.25_EB*(MU(I+1,J,K)+MU(I,J,K)+MU(I,J,K+1)+MU(I+1,J,K+1))
         TXZ(I,J,K) = MUY*(DUDZ + DWDX)
      ENDDO
   ENDDO
ENDDO
 
! Correct vorticity and stress tensor components at solid edges
 
EDGE_LOOP: DO IE=1,N_EDGES
   II  = IJKE(1,IE)
   JJ  = IJKE(2,IE)
   KK  = IJKE(3,IE)
   IEC = IJKE(4,IE)
   SELECT CASE(IEC)
      CASE(2)
         OMY(II,JJ,KK) = OME_E(IE)
         TXZ(II,JJ,KK) = TAU_E(IE)
   END SELECT
ENDDO EDGE_LOOP
 
! Compute gravity components
 
GX  = 0._EB
GZ  = EVALUATE_RAMP(T,DUMMY,I_RAMP_GZ)*GVEC(3)
 
! Upwind/Downwind bias factors
 
IF (PREDICTOR) THEN
   PMDT =  0.5_EB*DT
   MPDT = -0.5_EB*DT
ELSE
   PMDT = -0.5_EB*DT
   MPDT =  0.5_EB*DT
ENDIF
 
! Compute r-direction flux term FVX
 
IF (XS==0) THEN 
   I0 = 1
ELSE
   I0 = 0
ENDIF
 
J = 1
DO K= 1,KBAR
   DO I=I0,IBAR
      WP    = WW(I,J,K)   + WW(I+1,J,K)
      WM    = WW(I,J,K-1) + WW(I+1,J,K-1)
      EPSWP = 1._EB + WP*MPDT*RDZN(K)
      EPSWM = 1._EB + WM*PMDT*RDZN(K-1)
      WOMY  = EPSWP*WP*OMY(I,J,K) + EPSWM*WM*OMY(I,J,K-1)
      RRHO  = 2._EB/(RHOP(I,J,K)+RHOP(I+1,J,K))
      AH    = RHO_0(K)*RRHO - 1._EB   
      DWDZ  = (WW(I+1,J,K)-WW(I+1,J,K-1))*RDZ(K)
      TXXP  = MU(I+1,J,K)*( FOTH*DP(I+1,J,K) - 2._EB*DWDZ )
      DWDZ  = (WW(I,J,K)-WW(I,J,K-1))*RDZ(K)
      TXXM  = MU(I,J,K)  *( FOTH*DP(I,J,K) -2._EB*DWDZ )
      DTXXDX= RDXN(I)*(TXXP      -TXXM)
      DTXZDZ= RDZ(K) *(TXZ(I,J,K)-TXZ(I,J,K-1))
      DMUDX = (MU(I+1,J,K)-MU(I,J,K))*RDXN(I)
      VTRM  = RRHO*( DTXXDX + DTXZDZ - 2._EB*UU(I,J,K)*DMUDX/R(I) ) 
      FVX(I,J,K) = 0.25_EB*WOMY + GX*AH - VTRM 
   ENDDO
ENDDO   
 
! Compute z-direction flux term FVZ
 
J = 1
DO K=0,KBAR
   DO I=1,IBAR
      UP    = UU(I,J,K)   + UU(I,J,K+1)
      UM    = UU(I-1,J,K) + UU(I-1,J,K+1)
      EPSUP = 1._EB + UP*MPDT*RDXN(I)
      EPSUM = 1._EB + UM*PMDT*RDXN(I-1)
      UOMY  = EPSUP*UP*OMY(I,J,K) + EPSUM*UM*OMY(I-1,J,K)
      RRHO  = 2._EB/(RHOP(I,J,K)+RHOP(I,J,K+1))
      AH    = 0.5_EB*(RHO_0(K)+RHO_0(K+1))*RRHO - 1._EB
      DUDX  = (R(I)*UU(I,J,K+1)-R(I-1)*UU(I-1,J,K+1))*RDX(I)*RRN(I)
      TZZP  = MU(I,J,K+1)*( FOTH*DP(I,J,K+1) - 2._EB*DUDX )
      DUDX  = (R(I)*UU(I,J,K)-R(I-1)*UU(I-1,J,K))*RDX(I)*RRN(I)
      TZZM  = MU(I,J,K)  *( FOTH*DP(I,J,K)   - 2._EB*DUDX )
      DTXZDX= RDX(I) *(R(I)*TXZ(I,J,K)-R(I-1)*TXZ(I-1,J,K))*RRN(I)
      DTZZDZ= RDZN(K)*(TZZP      -TZZM)
      VTRM  = RRHO*(DTXZDX + DTZZDZ)
      FVZ(I,J,K) = -0.25_EB*UOMY + GZ*AH - VTRM 
   ENDDO
ENDDO   
 
! Baroclinic torque correction terms
 
IF (BAROCLINIC) CALL BAROCLINIC_CORRECTION
 
! Adjust FVX and FVZ at solid, internal obstructions for no flux
 
CALL NO_FLUX
 
END SUBROUTINE VELOCITY_FLUX_CYLINDRICAL
 
 
SUBROUTINE NO_FLUX

! Set FVX,FVY,FVZ inside internal blockages to maintain no flux

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP 
REAL(EB) :: RFODT
INTEGER  :: IC2,IC1,N,I,J,K,IW,II,JJ,KK,IOR
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
 
IF (PREDICTOR) THEN
   UU => U
   VV => V
   WW => W
ELSE
   UU => US
   VV => VS
   WW => WS
ENDIF
 
RFODT = RELAXATION_FACTOR/DT
 
! Drive velocity components inside solid obstructions towards zero
 
OBST_LOOP: DO N=1,N_OBST
   OB=>OBSTRUCTION(N)
   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         LOOP1: DO I=OB%I1  ,OB%I2
            IC1 = CELL_INDEX(I,J,K)
            IC2 = CELL_INDEX(I+1,J,K)
            IF (SOLID(IC1) .AND. SOLID(IC2)) FVX(I,J,K) = -RDXN(I)*(H(I+1,J,K)-H(I,J,K)) + RFODT*UU(I,J,K)
         ENDDO LOOP1
      ENDDO 
   ENDDO 
   DO K=OB%K1+1,OB%K2
      DO J=OB%J1  ,OB%J2
         LOOP2: DO I=OB%I1+1,OB%I2
            IC1 = CELL_INDEX(I,J,K)
            IC2 = CELL_INDEX(I,J+1,K)
            IF (SOLID(IC1) .AND. SOLID(IC2)) FVY(I,J,K) = -RDYN(J)*(H(I,J+1,K)-H(I,J,K)) + RFODT*VV(I,J,K)
         ENDDO LOOP2
      ENDDO 
   ENDDO 
   DO K=OB%K1  ,OB%K2
      DO J=OB%J1+1,OB%J2
         LOOP3: DO I=OB%I1+1,OB%I2
            IC1 = CELL_INDEX(I,J,K)
            IC2 = CELL_INDEX(I,J,K+1)
            IF (SOLID(IC1) .AND. SOLID(IC2)) FVZ(I,J,K) = -RDZN(K)*(H(I,J,K+1)-H(I,J,K)) + RFODT*WW(I,J,K)
         ENDDO LOOP3
      ENDDO 
   ENDDO 
ENDDO OBST_LOOP
 
! Add normal velocity to FVX, etc. for surface cells
 
WALL_LOOP: DO IW=1,NWC
 
!!! IF (.NOT.OBSTRUCTION(OBST_INDEX_W(IW))%SAWTOOTH) CYCLE WALL_LOOP  ! Allow velocity flux through surface of no SAWTOOTH obsts.

IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) CYCLE WALL_LOOP
IF (BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) CYCLE WALL_LOOP
II  = IJKW(1,IW)
JJ  = IJKW(2,IW)
KK  = IJKW(3,IW)
IOR = IJKW(4,IW)
SELECT CASE (BOUNDARY_TYPE(IW))
   CASE (SOLID_BOUNDARY,POROUS_BOUNDARY)
      SELECT CASE(IOR)
         CASE( 1) 
            FVX(II,JJ,KK)   = -RDXN(II)  *(H(II+1,JJ,KK)-H(II,JJ,KK)) + RFODT*(UU(II,JJ,KK)   + UWS(IW))
         CASE(-1) 
            FVX(II-1,JJ,KK) = -RDXN(II-1)*(H(II,JJ,KK)-H(II-1,JJ,KK)) + RFODT*(UU(II-1,JJ,KK) - UWS(IW))
         CASE( 2) 
            FVY(II,JJ,KK)   = -RDYN(JJ)  *(H(II,JJ+1,KK)-H(II,JJ,KK)) + RFODT*(VV(II,JJ,KK)   + UWS(IW))
         CASE(-2)
            FVY(II,JJ-1,KK) = -RDYN(JJ-1)*(H(II,JJ,KK)-H(II,JJ-1,KK)) + RFODT*(VV(II,JJ-1,KK) - UWS(IW))
         CASE( 3) 
            FVZ(II,JJ,KK)   = -RDZN(KK)  *(H(II,JJ,KK+1)-H(II,JJ,KK)) + RFODT*(WW(II,JJ,KK)   + UWS(IW))
         CASE(-3) 
            FVZ(II,JJ,KK-1) = -RDZN(KK-1)*(H(II,JJ,KK)-H(II,JJ,KK-1)) + RFODT*(WW(II,JJ,KK-1) - UWS(IW))
      END SELECT
   CASE (MIRROR_BOUNDARY)
      SELECT CASE(IOR)
         CASE( 1)
            FVX(II  ,JJ,KK) = 0._EB
         CASE(-1)
            FVX(II-1,JJ,KK) = 0._EB
         CASE( 2)
            FVY(II  ,JJ,KK) = 0._EB
         CASE(-2)
            FVY(II,JJ-1,KK) = 0._EB
         CASE( 3)
            FVZ(II  ,JJ,KK) = 0._EB
         CASE(-3)
            FVZ(II,JJ,KK-1) = 0._EB
      END SELECT
END SELECT 
 
ENDDO WALL_LOOP
 
END SUBROUTINE NO_FLUX
 
 
SUBROUTINE VELOCITY_PREDICTOR(T,NM,STOP_STATUS)
USE COMP_FUNCTIONS, ONLY: SECOND 
! Estimates the velocity components at the next time step
 
REAL(EB), INTENT(IN) :: T
REAL(EB) :: TNOW,RHS
INTEGER  :: STOP_STATUS,I,J,K
INTEGER, INTENT(IN) :: NM

IF (SOLID_PHASE_ONLY) RETURN
 
TNOW=SECOND() 
CALL POINT_TO_MESH(NM)
DO K=1,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         RHS = -FVX(I,J,K) - RDXN(I)*(H(I+1,J,K)-H(I,J,K))
         US(I,J,K) = U(I,J,K) + DT*RHS
      ENDDO 
   ENDDO 
ENDDO 

DO K=1,KBAR
   DO J=0,JBAR
      DO I=1,IBAR
         RHS = -FVY(I,J,K) - RDYN(J)*(H(I,J+1,K)-H(I,J,K))
         VS(I,J,K) = V(I,J,K) + DT*RHS
      ENDDO 
   ENDDO 
ENDDO 
 
DO K=0,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         RHS = -FVZ(I,J,K) - RDZN(K)*(H(I,J,K+1)-H(I,J,K))
         WS(I,J,K) = W(I,J,K) + DT*RHS
      ENDDO 
   ENDDO 
ENDDO 
! Check the stability criteria
 
DTOLD = DT
CALL CHECK_STABILITY(NM)
 
IF (DT < DTINT*1.E-4) THEN  ! The time step has gotten too small. Kill the job.
   STOP_STATUS = INSTABILITY_STOP
   RETURN
ENDIF
 
IF (.NOT.CHANGE_TIME_STEP(NM)) CALL VELOCITY_BC(T,NM)
 
TUSED(4,NM)=TUSED(4,NM)+SECOND()-TNOW
END SUBROUTINE VELOCITY_PREDICTOR
 
 
SUBROUTINE VELOCITY_CORRECTOR(T,NM)
USE COMP_FUNCTIONS, ONLY: SECOND
! Correct the velocity components
REAL(EB), INTENT(IN) :: T 
REAL(EB) :: TNOW,RHS
INTEGER  :: I,J,K
INTEGER, INTENT(IN) :: NM
 
IF (SOLID_PHASE_ONLY) RETURN

TNOW=SECOND() 
CALL POINT_TO_MESH(NM)

DO K=1,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         RHS = -FVX(I,J,K) - RDXN(I)*(H(I+1,J,K)-H(I,J,K))
         U(I,J,K) = .5_EB*(U(I,J,K) + US(I,J,K) + DT*RHS)
      ENDDO 
   ENDDO 
ENDDO 
DO K=1,KBAR
   DO J=0,JBAR
      DO I=1,IBAR
         RHS = -FVY(I,J,K) - RDYN(J)*(H(I,J+1,K)-H(I,J,K))
         V(I,J,K) = .5_EB*(V(I,J,K) + VS(I,J,K) + DT*RHS)
      ENDDO 
   ENDDO 
ENDDO 
 
DO K=0,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         RHS = -FVZ(I,J,K) - RDZN(K)*(H(I,J,K+1)-H(I,J,K))
         W(I,J,K) = .5_EB*(W(I,J,K) + WS(I,J,K) + DT*RHS)
      ENDDO 
   ENDDO 
ENDDO 

CALL VELOCITY_BC(T,NM)
 
TUSED(4,NM)=TUSED(4,NM)+SECOND()-TNOW
END SUBROUTINE VELOCITY_CORRECTOR
 
 
SUBROUTINE VELOCITY_BC(T,NM)

! Assert tangential velocity boundary conditions

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP 
REAL(EB) :: BC,MUA,T,FVT,UP,UM,VP,VM,WP,WM,DUDY,DUDZ,DVDX,DVDZ,DWDX,DWDY,TSI
INTEGER  :: I,J,K,IBC,NOM1,NOM2,IIO1,IIO2,JJO1,JJO2,KKO1,KKO2,NM,IE,II,JJ,KK,IEC,IOR,IWM,IWP,ICMM,ICMP,ICPM,ICPP,IC
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,U_Y,U_Z,V_X,V_Z,W_X,W_Y
TYPE (SURFACE_TYPE), POINTER :: SF

! Point to the appropriate velocity field

IF (PREDICTOR) THEN
   UU => US
   VV => VS
   WW => WS
ELSE
   UU => U
   VV => V
   WW => W
ENDIF

! Set the boundary velocity place holder to some large negative number

IF (CORRECTOR) THEN
   UVW_GHOST = -1.E6_EB
   U_Y => WORK1
   U_Z => WORK2
   V_X => WORK3
   V_Z => WORK4
   W_X => WORK5
   W_Y => WORK6
   U_Y = -1.E6_EB
   U_Z = -1.E6_EB
   V_X = -1.E6_EB
   V_Z = -1.E6_EB
   W_X = -1.E6_EB
   W_Y = -1.E6_EB
ENDIF

! Loop over all cell edges and determine the appropriate velocity BCs

EDGE_LOOP: DO IE=1,N_EDGES

   II   = IJKE( 1,IE)
   JJ   = IJKE( 2,IE)
   KK   = IJKE( 3,IE)
   IEC  = IJKE( 4,IE)
   IBC  = IJKE( 5,IE)
   IOR  = IJKE( 6,IE)
   NOM1 = IJKE( 7,IE)
   IIO1 = IJKE( 8,IE)
   JJO1 = IJKE( 9,IE)
   KKO1 = IJKE(10,IE)
   NOM2 = IJKE(11,IE)