part.f90

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

            CALL RANDOM_NUMBER(RN)
            THETA_RN = PY%SPRAY_ANGLE(1) + RN*(PY%SPRAY_ANGLE(2)-PY%SPRAY_ANGLE(1))
            CALL RANDOM_NUMBER(RN)
            PHI_RN = RN*TWOPI
         ENDIF PICK_PATTERN
         PHI_RN = PHI_RN + DV%ROTATION  ! Adjust for rotation of head by rotating about z-axis

         !  Adjust for tilt of sprinkler pipe

         SPHI   = SIN(PHI_RN)
         CPHI   = COS(PHI_RN)
         STHETA = SIN(THETA_RN)
         CTHETA = COS(THETA_RN)
         XTMP   = STHETA*CPHI
         YTMP   = STHETA*SPHI
         ZTMP   = -CTHETA
 
         ! First rotate about y-axis away from x-axis
 
         VLEN   = SQRT(XTMP**2+ZTMP**2)
         SHIFT1 = ACOS(ABS(XTMP)/VLEN)
            SELECT CASE(INT(SIGN(1._EB,ZTMP)))
            CASE (-1)
               IF (XTMP<0) SHIFT1 = PI-SHIFT1
            CASE ( 1)
            SELECT CASE(INT(SIGN(1._EB,XTMP)))
               CASE (-1)
                  SHIFT1 = SHIFT1+PI
               CASE ( 1)
                  SHIFT1 = TWOPI - SHIFT1
            END SELECT
         END SELECT
    
         SHIFT1 = SHIFT1 + TRIGT1
         XTMP = VLEN * COS(SHIFT1)
         ZTMP = -VLEN * SIN(SHIFT1)
 
         ! Second rotate about z-axis away from x-axis
 
         VLEN   = SQRT(XTMP**2+YTMP**2)
         SHIFT1 = ACOS(ABS(XTMP)/VLEN)
         SELECT CASE(INT(SIGN(1._EB,YTMP)))
            CASE ( 1)
               IF (XTMP<0) SHIFT1 = PI-SHIFT1
            CASE (-1)
            SELECT CASE(INT(SIGN(1._EB,XTMP)))
               CASE (-1)
                  SHIFT1 = SHIFT1+PI
               CASE ( 1) 
                  SHIFT1 = TWOPI - SHIFT1
            END SELECT
         END SELECT
 
         SHIFT2 = TRIGT2
         SELECT CASE(INT(SIGN(1._EB,DV%ORIENTATION(1))))
            CASE (-1)
               IF (DV%ORIENTATION(2)>0) SHIFT2 = TWOPI - SHIFT2
            CASE ( 1)
            SELECT CASE(INT(SIGN(1._EB,DV%ORIENTATION(2))))
               CASE (-1) 
                  SHIFT2 = PI-SHIFT2
               CASE ( 1)
                  SHIFT2 = SHIFT2+ PI
            END SELECT
         END SELECT
         SHIFT1=SHIFT1+SHIFT2
         XTMP = VLEN * COS(SHIFT1)
         YTMP = VLEN * SIN(SHIFT1)
         DROPLET_SPEED = PY%DROPLET_VELOCITY
 
         ! Compute initial position and velocity of droplets
 
         DR%U = DROPLET_SPEED*XTMP
         DR%V = DROPLET_SPEED*YTMP
         DR%W = DROPLET_SPEED*ZTMP
         DR%X = DV%X + PY%OFFSET*XTMP
         DR%Y = DV%Y + PY%OFFSET*YTMP
         DR%Z = DV%Z + PY%OFFSET*ZTMP
         IF (TWO_D) THEN
            DR%V = 0._EB
            DR%Y = DV%Y
         ENDIF
         IF (DR%X<=XS .OR. DR%X>=XF) CYCLE CHOOSE_COORDS
         IF (DR%Y<=YS .OR. DR%Y>=YF) CYCLE CHOOSE_COORDS
         IF (DR%Z<=ZS .OR. DR%Z>=ZF) CYCLE CHOOSE_COORDS
         XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
         YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
         ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
         II = FLOOR(XI+1._EB)
         JJ = FLOOR(YJ+1._EB)
         KK = FLOOR(ZK+1._EB)
         IC = CELL_INDEX(II,JJ,KK)
         IF (.NOT.SOLID(IC)) EXIT CHOOSE_COORDS
 
      ENDDO CHOOSE_COORDS

      ! Randomly choose droplet size according to Cumulative Distribution Function (CDF)
 
      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
 
      ! Sum up mass of liquid being introduced

      IF(PY%SPRAY_PATTERN_INDEX>0) THEN
         MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3*SOLID_ANGLE_FLOWRATE
      ELSE
         MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3
      ENDIF
      DROP_SUM = DROP_SUM + 1
   ENDDO DROPLET_INSERT_LOOP
 
   ! Compute weighting factor for the droplets just inserted

   IF (DROP_SUM > 0) THEN
      IF(PY%SPRAY_PATTERN_INDEX>0) THEN
         PWT0 = FLOW_RATE*PC%DT_INSERT/(MASS_SUM*REAL(PC%N_INSERT,EB)/SOLID_ANGLE_TOTAL_FLOWRATE)
      ELSE
         PWT0 = FLOW_RATE*PC%DT_INSERT/MASS_SUM
      ENDIF
      DROPLET(NLP-DROP_SUM+1:NLP)%PWT = DROPLET(NLP-DROP_SUM+1:NLP)%PWT*PWT0
   ENDIF

   ! Indicate that droplets from this device have been inserted at this time T

   DV%T = T 

ENDDO SPRINKLER_INSERT_LOOP
 
! Loop through all boundary cells and insert particles if appropriate
 
WALL_INSERT_LOOP: DO IW=1,NWC
 
   IF (T < TW(IW))                       CYCLE WALL_INSERT_LOOP
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) CYCLE WALL_INSERT_LOOP
   IBC =  IJKW(5,IW)
   SF  => SURFACE(IBC)
   IPC =  SF%PART_INDEX
   IF (IPC < 1)    CYCLE WALL_INSERT_LOOP
   PC  => PARTICLE_CLASS(IPC)
   IF (T < PC%INSERT_CLOCK(NM)) CYCLE WALL_INSERT_LOOP
   IF (UW(IW) >= -0.0001_EB) CYCLE WALL_INSERT_LOOP
 
   II = IJKW(1,IW)
   JJ = IJKW(2,IW)
   KK = IJKW(3,IW)
   IC = CELL_INDEX(II,JJ,KK)
   IF (.NOT.SOLID(IC)) CYCLE WALL_INSERT_LOOP
 
   IF (NM > 1) THEN
      IIG = IJKW(6,IW)
      JJG = IJKW(7,IW)
      KKG = IJKW(8,IW)
      IF (INTERPOLATED_MESH(IIG,JJG,KKG) > 0) CYCLE WALL_INSERT_LOOP
   ENDIF
 
   IOR = IJKW(4,IW)
   MASS_SUM = 0._EB

   PARTICLE_INSERT_LOOP: DO I=1,NPPCW(IW)  ! Loop over all particles for the IW'th cell
 
      IF (NLP+1 > MAXIMUM_DROPLETS) EXIT PARTICLE_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)
 
      IF (ABS(IOR)==1) THEN
         IF (IOR== 1) DR%X = X(II)   + VENT_OFFSET*DX(II+1)
         IF (IOR==-1) DR%X = X(II-1) - VENT_OFFSET*DX(II-1)
         CALL RANDOM_NUMBER(RN)
         DR%Y = Y(JJ-1) + DY(JJ)*RN
         CALL RANDOM_NUMBER(RN)
         DR%Z = Z(KK-1) + DZ(KK)*RN
      ENDIF
      IF (ABS(IOR)==2) THEN
         IF (IOR== 2) DR%Y = Y(JJ)   + VENT_OFFSET*DY(JJ+1)
         IF (IOR==-2) DR%Y = Y(JJ-1) - VENT_OFFSET*DY(JJ-1)
         CALL RANDOM_NUMBER(RN)
         DR%X = X(II-1) + DX(II)*RN
         CALL RANDOM_NUMBER(RN)
         DR%Z = Z(KK-1) + DZ(KK)*RN
      ENDIF
      IF (ABS(IOR)==3) THEN
         IF (IOR== 3) DR%Z = Z(KK)   + VENT_OFFSET*DZ(KK+1)
         IF (IOR==-3) DR%Z = Z(KK-1) - VENT_OFFSET*DZ(KK-1)
         CALL RANDOM_NUMBER(RN)
         DR%X = X(II-1) + DX(II)*RN
         CALL RANDOM_NUMBER(RN)
         DR%Y = Y(JJ-1) + DY(JJ)*RN
      ENDIF
 
      SELECT CASE(IOR)  ! Give particles an initial velocity (if desired)
         CASE( 1) 
            DR%U = -UW(IW)
         CASE(-1) 
            DR%U =  UW(IW)
         CASE( 2) 
            DR%V = -UW(IW)
         CASE(-2) 
            DR%V =  UW(IW)
         CASE( 3)
            DR%W = -UW(IW)
         CASE(-3) 
            DR%W =  UW(IW)
      END SELECT
 
      ! Save the insertion time (TP) and scalar property (SP) for the particle
 
      DR%A_X    = 0._EB
      DR%A_Y    = 0._EB
      DR%A_Z    = 0._EB
      DR%CLASS  = IPC
      DR%TAG    = PARTICLE_TAG
      IF (MOD(NLP,PC%SAMPLING)==0) THEN
         DR%SHOW = .TRUE.
      ELSE
         DR%SHOW = .FALSE.
      ENDIF
      DR%R   = 0.0_EB
      DR%TMP = PC%TMP_INITIAL
      DR%T   = T
      DR%PWT = 1._EB
      DR%IOR = 0
 
      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 PARTICLE_INSERT_LOOP

   IF (MASS_SUM > 0._EB) THEN
      IF (SF%PARTICLE_MASS_FLUX > 0._EB) THEN
         DROPLET(NLP-NPPCW(IW)+1:NLP)%PWT = &
         DROPLET(NLP-NPPCW(IW)+1:NLP)%PWT*SF%PARTICLE_MASS_FLUX*AREA_ADJUST(IW)*AW(IW)*PC%DT_INSERT/MASS_SUM
      ENDIF
   ENDIF

ENDDO WALL_INSERT_LOOP

! Reset particle/droplet insertion clock
 
DO IPC=1,N_PART
   PC => PARTICLE_CLASS(IPC)
   IF (T >= PC%INSERT_CLOCK(NM)) PC%INSERT_CLOCK(NM) = PC%INSERT_CLOCK(NM) + PC%DT_INSERT
   IF (T >= PC%INSERT_CLOCK(NM)) INSERT_ANOTHER_BATCH = .TRUE.
ENDDO

IF (.NOT.INSERT_ANOTHER_BATCH) EXIT OVERALL_INSERT_LOOP
ENDDO OVERALL_INSERT_LOOP
 
TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE INSERT_DROPLETS_AND_PARTICLES
 
 

SUBROUTINE RANDOM_CHOICE(CDF,VAR,NPTS,CHOICE)
 
INTEGER,  INTENT(IN)  :: NPTS
REAL(EB), INTENT(IN)  :: CDF(0:NPTS),VAR(0:NPTS)
REAL(EB), INTENT(OUT) :: CHOICE
INTEGER  :: IT
REAL(EB) :: CFRAC,RN,A,B
 
CALL RANDOM_NUMBER(RN)
A = MINVAL(CDF)
B = MAXVAL(CDF)
RN = A + (B-A)*RN
 
CDF_LOOP: DO IT=1,NPTS
   IF (CDF(IT) > RN) THEN
      CFRAC  = (RN-CDF(IT-1))/(CDF(IT)-CDF(IT-1))
      CHOICE = VAR(IT-1) + (VAR(IT)-VAR(IT-1))*CFRAC
      EXIT CDF_LOOP
   ENDIF
ENDDO CDF_LOOP
 
END SUBROUTINE RANDOM_CHOICE
 
 

SUBROUTINE UPDATE_PARTICLES(T,NM)

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

IF (MESHES(NM)%NLP==0) RETURN
IF (EVACUATION_ONLY(NM)) RETURN
TNOW=SECOND()
CALL POINT_TO_MESH(NM)

IF (CORRECTOR) CALL MOVE_PARTICLES(T,NM)
IF (CORRECTOR) CALL PARTICLE_MASS_ENERGY_TRANSFER(T,NM)
CALL PARTICLE_MOMENTUM_TRANSFER

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



SUBROUTINE MOVE_PARTICLES(T,NM)

! Momentum and heat transfer from all particles and droplets
 
USE COMP_FUNCTIONS, ONLY : SECOND  
USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP, AFILL2
USE PHYSICAL_FUNCTIONS, ONLY : DRAG
 
REAL(EB) :: RHO_G,RVC,RDS,RDC,QREL,SFAC,UREL,VREL,WREL,TMP_G,RN,THETA_RN, &
            RD,T,C_DRAG,XI,YJ,ZK,MU_AIR, &
            DTSP,DTMIN,UBAR,VBAR,WBAR,BFAC,GRVT1,GRVT2,GRVT3,AUREL,AVREL,AWREL,CONST, &
            UVW,DUMMY=0._EB,X_OLD,Y_OLD,Z_OLD,STEP_FRACTION(-3:3),R_NEW,SURFACE_DROPLET_DIAMETER
LOGICAL :: HIT_SOLID
INTEGER :: ICN,I,IIN,JJN,KKN,II,JJ,KK,IIX,JJY,KKZ,IW,N,NITER,IWP1,IWM1,IWP2,IWM2,IWP3,IWM3,IOR_OLD,IC,IOR_FIRST,IML
INTEGER, INTENT(IN) :: NM

SURFACE_DROPLET_DIAMETER = 0.001_EB  ! All droplets adjusted to this size when on solid (m)

GRVT1 = -EVALUATE_RAMP(T,DUMMY,I_RAMP_GX)*GVEC(1) 
GRVT2 = -EVALUATE_RAMP(T,DUMMY,I_RAMP_GY)*GVEC(2) 
GRVT3 = -EVALUATE_RAMP(T,DUMMY,I_RAMP_GZ)*GVEC(3) 

! Loop through all Lagrangian particles and move them one time step

DROPLET_LOOP: DO I=1,NLP  

   DR => DROPLET(I)
   PC => PARTICLE_CLASS(DR%CLASS)

   ! Determine the current coordinates of the particle

   XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
   YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
   ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
   II  = FLOOR(XI+1._EB)
   JJ  = FLOOR(YJ+1._EB)
   KK  = FLOOR(ZK+1._EB)
   IC  = CELL_INDEX(II,JJ,KK)

   ! Throw out particles that are inside a solid obstruction

   IF (SOLID(IC)) THEN
      DR%X = 1.E6_EB
      CYCLE DROPLET_LOOP
   ENDIF

   ! Interpolate the nearest velocity components of the gas

   IIX  = FLOOR(XI+.5_EB)
   JJY  = FLOOR(YJ+.5_EB)
   KKZ  = FLOOR(ZK+.5_EB)
   UBAR = AFILL2(U,II-1,JJY,KKZ,XI-II+1,YJ-JJY+.5_EB,ZK-KKZ+.5_EB)
   VBAR = AFILL2(V,IIX,JJ-1,KKZ,XI-IIX+.5_EB,YJ-JJ+1,ZK-KKZ+.5_EB)
   WBAR = AFILL2(W,IIX,JJY,KK-1,XI-IIX+.5_EB,YJ-JJY+.5_EB,ZK-KK+1)
    
   ! If the particle is just a massless tracer, just move it and go on to the next particle

   IF (PC%MASSLESS) THEN
      DR%U = UBAR
      DR%V = VBAR
      DR%W = WBAR
      DR%X = DR%X + DR%U*DT
      DR%Y = DR%Y + DR%V*DT
      DR%Z = DR%Z + DR%W*DT
      CYCLE DROPLET_LOOP
   ENDIF

   ! Calculate the particle velocity components and the amount of momentum to transfer to the gas

   RVC = RDX(II)*RDY(JJ)*RDZ(KK)
   RD  = DR%R
   RDS = RD*RD
   RDC = RD*RDS

   UREL  = DR%U - UBAR