FDS Source: add experimental model for ember snag (sticking to other …

…particles), activated with EMBER_SNAG_FACTOR

Source/part.f90

@@ -2520,11 +2520,14 @@ SUBROUTINE MOVE_IN_GAS
             DD,DD_X,DD_Y,DD_Z,DW_X,DW_Y,DW_Z,K_TERM(3),Y_TERM(3),C_DRAG,A_DRAG,&
             GX_LOC,GY_LOC,GZ_LOC,DRAG_MAX(3)=0._EB,K_SGS,U_P,KN,M_DOT,&
             EMBER_DENSITY,EMBER_VOLUME,ACCEL_X,ACCEL_Y,ACCEL_Z,&
-            LP_FORCE,FACE_VOLS(2,2,2),VEL_G_INT(3),VOL_WGT(2,2,2)
+            LP_FORCE,FACE_VOLS(2,2,2),VEL_G_INT(3),VOL_WGT(2,2,2),&
+            EMBER_PACKING_RATIO,LOCAL_PACKING_RATIO,LPC_GEOM_FACTOR
 REAL(EB) :: WGT(2,2,2,3)
 REAL(EB), POINTER, DIMENSION(:,:,:) :: FV_D=>NULL(),VEL_G=>NULL()
 REAL(EB), SAVE :: BETA
-INTEGER :: IIX,JJY,KKZ,IL,JL,KL,AXIS
+INTEGER :: IIX,JJY,KKZ,IL,JL,KL,AXIS,N_LPC2
+LOGICAL :: STUCK
+TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC2
 
 ! Save current values of particle velocity components
 
@@ -2640,6 +2643,49 @@ SUBROUTINE MOVE_IN_GAS
    WBAR = WBAR + U_P*DW_Z
 ENDIF
 
+IF (LPC%EMBER_PARTICLE) THEN
+   SELECT CASE(SF%GEOMETRY)
+      CASE(SURF_CARTESIAN)
+         EMBER_VOLUME = SF%LENGTH * SF%WIDTH * 2._EB*R_D
+      CASE(SURF_CYLINDRICAL)
+         EMBER_VOLUME = SF%LENGTH * PI*R_D**2
+      CASE(SURF_SPHERICAL)
+         EMBER_VOLUME = FOTHPI * R_D**3
+   END SELECT
+
+   ! experimental ember snag model
+   IF (LP%EMBER .AND. LPC%EMBER_SNAG_FACTOR>0._EB) THEN
+      STUCK=.FALSE.
+      ! constrain max packing ratio to 1
+      EMBER_PACKING_RATIO = MIN(1._EB, EMBER_VOLUME * LP%PWT * LP%RVC)
+      LOCAL_PACKING_RATIO = 0._EB
+      LPC2_LOOP: DO N_LPC2=1,N_LAGRANGIAN_CLASSES
+         LPC2 => LAGRANGIAN_PARTICLE_CLASS(N_LPC2)
+         IF (LPC2%ARRAY_INDEX == LPC%ARRAY_INDEX) CYCLE LPC2_LOOP
+         IF (AVG_DROP_RAD(IIG_OLD,JJG_OLD,KKG_OLD,LPC2%ARRAY_INDEX)<TWO_EPSILON_EB) CYCLE LPC2_LOOP
+         SELECT CASE(SURFACE(LPC2%SURF_INDEX)%GEOMETRY)
+            CASE(SURF_CARTESIAN)
+               LPC_GEOM_FACTOR = 1._EB
+            CASE(SURF_CYLINDRICAL)
+               LPC_GEOM_FACTOR = 0.5_EB*PI
+            CASE(SURF_SPHERICAL)
+               LPC_GEOM_FACTOR = FOTH
+         END SELECT
+         LOCAL_PACKING_RATIO = LOCAL_PACKING_RATIO + AVG_DROP_AREA(IIG_OLD,JJG_OLD,KKG_OLD,LPC2%ARRAY_INDEX)*&
+                  AVG_DROP_RAD(IIG_OLD,JJG_OLD,KKG_OLD,LPC2%ARRAY_INDEX)*LPC_GEOM_FACTOR
+      ENDDO LPC2_LOOP
+      LOCAL_PACKING_RATIO = MIN(1._EB, LOCAL_PACKING_RATIO)
+      CALL RANDOM_NUMBER(RN)
+      IF (LOCAL_PACKING_RATIO>TWO_EPSILON_EB .AND. &
+         RN<(LOCAL_PACKING_RATIO*EMBER_PACKING_RATIO)**LPC%EMBER_SNAG_FACTOR) THEN
+         STUCK=.TRUE.
+         LP%U = 0._EB
+         LP%V = 0._EB
+         LP%W = 0._EB
+      ENDIF
+   ENDIF
+ENDIF
+
 ! Calculate the particle drag coefficient
 
 RHO_G = RHO(IIG_OLD,JJG_OLD,KKG_OLD)
@@ -2760,14 +2806,6 @@ SUBROUTINE MOVE_IN_GAS
 ! Experimental ember generation model
 
 IF (LPC%EMBER_PARTICLE .AND. .NOT.LP%EMBER) THEN
-   SELECT CASE(SF%GEOMETRY)
-      CASE(SURF_CARTESIAN)
-         EMBER_VOLUME = SF%LENGTH * SF%WIDTH * 2._EB*R_D
-      CASE(SURF_CYLINDRICAL)
-         EMBER_VOLUME = SF%LENGTH * PI*R_D**2
-      CASE(SURF_SPHERICAL)
-         EMBER_VOLUME = FOTHPI * R_D**3
-   END SELECT
    EMBER_DENSITY = 0._EB
    IF (EMBER_VOLUME>TWO_EPSILON_EB) EMBER_DENSITY = LP%MASS/EMBER_VOLUME
    IF ( EMBER_DENSITY < LPC%EMBER_DENSITY_THRESHOLD .AND. &
@@ -2782,7 +2820,7 @@ SUBROUTINE MOVE_IN_GAS
 
 ! Move the particles unless they are STATIC
 
-PARTICLE_NON_STATIC_IF: IF (.NOT.LPC%STATIC .OR. LP%EMBER) THEN ! Move airborne, non-stationary particles
+PARTICLE_NON_STATIC_IF: IF (.NOT.LPC%STATIC .OR. (LP%EMBER .AND. .NOT.STUCK)) THEN ! Move airborne, non-stationary particles
 
    ! Compute gravity components
 

Source/read.f90

@@ -5572,7 +5572,7 @@ SUBROUTINE READ_PART
             SURFACE_TENSION,BREAKUP_RATIO,BREAKUP_GAMMA_D,BREAKUP_SIGMA_D,&
             DENSE_VOLUME_FRACTION,REAL_REFRACTIVE_INDEX,COMPLEX_REFRACTIVE_INDEX,RUNNING_AVERAGE_FACTOR,&
             RUNNING_AVERAGE_FACTOR_WALL,KILL_DIAMETER,&
-            EMBER_DENSITY_THRESHOLD,EMBER_VELOCITY_THRESHOLD,PRIMARY_BREAKUP_LENGTH,&
+            EMBER_DENSITY_THRESHOLD,EMBER_VELOCITY_THRESHOLD,EMBER_SNAG_FACTOR,PRIMARY_BREAKUP_LENGTH,&
             PRIMARY_BREAKUP_DRAG_REDUCTION_FACTOR,HEAT_TRANSFER_COEFFICIENT_GAS,HEAT_TRANSFER_COEFFICIENT_SOLID,&
             MASS_TRANSFER_COEFFICIENT
 REAL(EB) :: DRAG_COEFFICIENT(3),FREE_AREA_FRACTION,PERMEABILITY(3),POROUS_VOLUME_FRACTION,SHAPE_FACTOR
@@ -5587,7 +5587,7 @@ SUBROUTINE READ_PART
                 BREAKUP_SIGMA_D,CHECK_DISTRIBUTION,CNF_RAMP_ID,COLOR,COMPLEX_REFRACTIVE_INDEX,&
                 CTRL_ID,DEBUG,DENSE_VOLUME_FRACTION,&
                 DEVC_ID,DIAMETER,DISTRIBUTION,DRAG_COEFFICIENT,DRAG_LAW,&
-                EMBER_DENSITY_THRESHOLD,EMBER_PARTICLE,EMBER_VELOCITY_THRESHOLD,EVAP_MODEL,&
+                EMBER_DENSITY_THRESHOLD,EMBER_PARTICLE,EMBER_VELOCITY_THRESHOLD,EMBER_SNAG_FACTOR,EVAP_MODEL,&
                 FREE_AREA_FRACTION,FYI,GAMMA_D,HEAT_OF_COMBUSTION,HEAT_TRANSFER_COEFFICIENT_GAS,HEAT_TRANSFER_COEFFICIENT_SOLID,&
                 HORIZONTAL_VELOCITY,ID,INITIAL_TEMPERATURE,KILL_DIAMETER,MASSLESS,&
                 MASS_TRANSFER_COEFFICIENT,MAXIMUM_DIAMETER,&
@@ -5917,6 +5917,7 @@ SUBROUTINE READ_PART
    LPC%TRACK_EMBERS        = TRACK_EMBERS
    LPC%EMBER_DENSITY_THRESHOLD = EMBER_DENSITY_THRESHOLD
    LPC%EMBER_VELOCITY_THRESHOLD = EMBER_VELOCITY_THRESHOLD
+   LPC%EMBER_SNAG_FACTOR   = EMBER_SNAG_FACTOR
    LPC%PRIMARY_BREAKUP_TIME = PRIMARY_BREAKUP_LENGTH ! user enters LENGTH, later divide by PARTICLE_VELOCITY to get TIME
    LPC%PRIMARY_BREAKUP_DRAG_REDUCTION_FACTOR = PRIMARY_BREAKUP_DRAG_REDUCTION_FACTOR
 
@@ -6114,6 +6115,7 @@ SUBROUTINE SET_PART_DEFAULTS
 TRACK_EMBERS             = .TRUE.
 EMBER_DENSITY_THRESHOLD  = 0._EB
 EMBER_VELOCITY_THRESHOLD = 1000._EB
+EMBER_SNAG_FACTOR        = -1._EB
 PRIMARY_BREAKUP_LENGTH   = -1._EB
 PRIMARY_BREAKUP_DRAG_REDUCTION_FACTOR = 1._EB
 POROUS_VOLUME_FRACTION = -1._EB
@@ -13130,7 +13132,7 @@ SUBROUTINE READ_DEVC
                             ' XYZ will be set to the center of XB.'
       IF (MY_RANK==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
    ENDIF
-   
+
    ! Check the QUANTITY_RANGE
 
    IF (QUANTITY_RANGE(2) <= QUANTITY_RANGE(1)) THEN
@@ -15153,7 +15155,7 @@ SUBROUTINE READ_SLCF
       IF (ABS(XB(1)-XB(2))<TWO_EPSILON_EB) IOR = 1
       IF (ABS(XB(3)-XB(4))<TWO_EPSILON_EB) IOR = 2
       IF (ABS(XB(5)-XB(6))<TWO_EPSILON_EB) IOR = 3
-      
+
       CALL CHECK_XB(XB)
 
       XB(1) = MAX(XB(1),XS)
@@ -15187,14 +15189,14 @@ SUBROUTINE READ_SLCF
             ENDIF
          CASE DEFAULT
             IF (XB(1)>XF .OR. XB(2)<XS .OR. XB(3)>YF .OR. XB(4)<YS .OR. XB(5)>ZF .OR. XB(6)<ZS) CULL_SLICE = .TRUE.
-      END SELECT      
+      END SELECT
       IF (CULL_SLICE) THEN
          N_SLCF = N_SLCF - 1
          IF (VECTOR .AND. TWO_D) N_SLCF = N_SLCF - 2
          IF (VECTOR .AND. .NOT. TWO_D) N_SLCF = N_SLCF - 3
          CYCLE SLCF_LOOP
       ENDIF
-      
+
       ! Process vector quantities
 
       NITER = 1

Source/type.f90

@@ -90,6 +90,7 @@ MODULE TYPES
    REAL(EB) :: SHAPE_FACTOR               !< Ratio of particle cross sectional area to surface area
    REAL(EB) :: EMBER_DENSITY_THRESHOLD    !< Density at which vegetative particle becomes a flying ember
    REAL(EB) :: EMBER_VELOCITY_THRESHOLD   !< Velocity at which vegetative particle becomes a flying ember
+   REAL(EB) :: EMBER_SNAG_FACTOR          !< Scaling factor for probability of ember snaging in a collection of particles
    REAL(EB) :: PRIMARY_BREAKUP_TIME       !< Time (s) after insertion when droplet breaks up
    REAL(EB) :: PRIMARY_BREAKUP_DRAG_REDUCTION_FACTOR   !< Drag reduction factor
    REAL(EB) :: RUNNING_AVERAGE_FACTOR_WALL             !< Fraction of old value used in summations of droplets stuck to walls