Skip to content

Commit

Permalink
Merge pull request #12225 from mcgratta/master
Browse files Browse the repository at this point in the history 
FDS Source: Remove SF_FIXED from PARTICLE_RUNNING_AVERAGES and do som…
  • Loading branch information
@mcgratta
mcgratta authored Nov 12, 2023
2 parents b90adf2 + 003127f commit d8822ca
Showing 1 changed file with 139 additions and 146 deletions.
285 changes: 139 additions & 146 deletions Source/part.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2962,14 +2962,6 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)
REAL(EB), INTENT(IN) :: T,DT
INTEGER, INTENT(IN) :: NM

REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_DEN=>NULL()
!< Average particle density in a grid cell (kg/m3) used in updating AVG_DROP_DEN
REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_RAD=>NULL()
!< Average particle radius in a grid cell (m) used in updating AVG_DROP_RAD
REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_TMP=>NULL()
!< Average particle temperature in a grid cell (K) used in updating AVG_DROP_TMP
REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_AREA=>NULL()
!< Average particle cross sectional areat in a grid cell (m2) used in updating AVG_DROP_AREA
REAL(EB), POINTER, DIMENSION(:,:,:) :: MVAP_TOT=>NULL()
!< Amount of mass evaporated into a grid cell (kg)
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_INTERIM=>NULL()
Expand Down Expand Up @@ -3031,8 +3023,6 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)
REAL(EB) :: SH_FAC_GAS !< Sherwood number of a particle in the gas
REAL(EB) :: M_VAP !< Mass evaporated (kg) from the particle in the current sub time step
REAL(EB) :: M_VAP_MAX !< Maximum allowable evaporation (kg)
REAL(EB) :: DEN_ADD !< Weighted drop mass divided by cell volume (kg/m3)
REAL(EB) :: AREA_ADD !< Weighted drop area divided by cell volume (1/m)
REAL(EB) :: Y_ALL(1:N_SPECIES) !< Mass fraction of all primitive species
REAL(EB) :: X_DROP !< Equilibirum vapor mole fraction at the particle temperature
REAL(EB) :: Y_DROP !< Equilibrium vapor mass fraction at the particle temperature
Expand Down Expand Up @@ -3112,7 +3102,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)
CALL POINT_TO_MESH(NM)

IF (MESHES(NM)%NLP==0) THEN
CALL PARTICLE_RUNNING_AVERAGES
IF (N_LP_ARRAY_INDICES > 0) CALL PARTICLE_RUNNING_AVERAGES
T_USED(8)=T_USED(8)+CURRENT_TIME()-TNOW
RETURN
ELSE
Expand Down Expand Up @@ -3805,7 +3795,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)

! Sum up various quantities used in running averages

CALL PARTICLE_RUNNING_AVERAGES
IF (N_LP_ARRAY_INDICES > 0) CALL PARTICLE_RUNNING_AVERAGES

! Remove PARTICLEs that have completely evaporated

Expand All @@ -3820,170 +3810,173 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)

SUBROUTINE PARTICLE_RUNNING_AVERAGES

SUM_PART_QUANTITIES: IF (N_LP_ARRAY_INDICES > 0) THEN
INTEGER :: I,J,K
REAL(EB) :: DEN_ADD,AREA_ADD
REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_DEN,DROP_RAD,DROP_TMP,DROP_AREA

DROP_AREA => WORK1
DROP_DEN => WORK4
DROP_RAD => WORK5
DROP_TMP => WORK6
DROP_AREA => WORK1
DROP_DEN => WORK4
DROP_RAD => WORK5
DROP_TMP => WORK6

PART_CLASS_SUM_LOOP: DO N_LPC = 1,N_LAGRANGIAN_CLASSES
PART_CLASS_SUM_LOOP: DO N_LPC = 1,N_LAGRANGIAN_CLASSES

LPC => LAGRANGIAN_PARTICLE_CLASS(N_LPC)
IF (LPC%MASSLESS_TRACER .OR. LPC%MASSLESS_TARGET) CYCLE PART_CLASS_SUM_LOOP
LPC => LAGRANGIAN_PARTICLE_CLASS(N_LPC)
IF (LPC%MASSLESS_TRACER .OR. LPC%MASSLESS_TARGET) CYCLE PART_CLASS_SUM_LOOP

DROP_DEN = 0._EB
DROP_TMP = 0._EB
DROP_RAD = 0._EB
DROP_AREA = 0._EB
DROP_DEN = 0._EB
DROP_TMP = 0._EB
DROP_RAD = 0._EB
DROP_AREA = 0._EB

PARTICLE_LOOP_2: DO IP=1,NLP
PARTICLE_LOOP_2: DO IP=1,NLP

LP => LAGRANGIAN_PARTICLE(IP)
IF (LP%CLASS_INDEX /= N_LPC) CYCLE PARTICLE_LOOP_2
BC => BOUNDARY_COORD(LP%BC_INDEX)
LP => LAGRANGIAN_PARTICLE(IP)
IF (LP%CLASS_INDEX /= N_LPC) CYCLE PARTICLE_LOOP_2
BC => BOUNDARY_COORD(LP%BC_INDEX)

! Determine the mass of the PARTICLE/particle, depending on whether the particle has a distinct SURFace type.

IF (LPC%LIQUID_DROPLET) THEN
R_DROP = LP%RADIUS
A_DROP = PI*R_DROP**2
ELSE
SF => SURFACE(LPC%SURF_INDEX)
R_DROP = LP%RADIUS
SELECT CASE(SF%GEOMETRY)
CASE(SURF_CARTESIAN)
A_DROP = 2._EB*SF%LENGTH*SF%WIDTH
CASE(SURF_CYLINDRICAL)
A_DROP = 2._EB*SF%LENGTH*R_DROP
CASE(SURF_SPHERICAL)
A_DROP = PI*R_DROP**2
END SELECT
ENDIF

! Assign particle or PARTICLE mass to the grid cell if the particle/PARTICLE not on a surface

IF (BC%IOR==0 .OR. LPC%SOLID_PARTICLE) THEN
LP_B1 => BOUNDARY_PROP1(LP%B1_INDEX)
II = BC%IIG
JJ = BC%JJG
KK = BC%KKG
DEN_ADD = LP%PWT*LP%MASS*LP%RVC
AREA_ADD = LP%PWT*A_DROP*LP%RVC
DROP_DEN(BC%IIG,BC%JJG,BC%KKG) = DROP_DEN(BC%IIG,BC%JJG,BC%KKG) + DEN_ADD
DROP_TMP(BC%IIG,BC%JJG,BC%KKG) = DROP_TMP(BC%IIG,BC%JJG,BC%KKG) + DEN_ADD*LP_B1%TMP_F
DROP_RAD(BC%IIG,BC%JJG,BC%KKG) = DROP_RAD(BC%IIG,BC%JJG,BC%KKG) + AREA_ADD*R_DROP
DROP_AREA(BC%IIG,BC%JJG,BC%KKG) = DROP_AREA(BC%IIG,BC%JJG,BC%KKG) + AREA_ADD
ENDIF

! Determine the mass of the PARTICLE/particle, depending on whether the particle has a distinct SURFace type.
! Compute Mass Per Unit Area (MPUA) for a liquid droplet stuck to a solid surface

IF (BC%IOR/=0 .AND. (LP%WALL_INDEX>0 .OR. LP%CFACE_INDEX>0)) THEN
IF (LP%WALL_INDEX>0) THEN
WC => WALL(LP%WALL_INDEX)
B1 => BOUNDARY_PROP1(WC%B1_INDEX)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
ELSE
CFA => CFACE(LP%CFACE_INDEX)
B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
ENDIF
IF (LPC%LIQUID_DROPLET) THEN
R_DROP = LP%RADIUS
A_DROP = PI*R_DROP**2
FTPR = FOTHPI * LPC%DENSITY
M_DROP = FTPR*R_DROP**3
ELSE
SF => SURFACE(LPC%SURF_INDEX)
R_DROP = LP%RADIUS
SELECT CASE(SF%GEOMETRY)
CASE(SURF_CARTESIAN)
A_DROP = 2._EB*SF%LENGTH*SF%WIDTH
CASE(SURF_CYLINDRICAL)
A_DROP = 2._EB*SF%LENGTH*R_DROP
CASE(SURF_SPHERICAL)
A_DROP = PI*R_DROP**2
END SELECT
M_DROP = LP%MASS
ENDIF
B2%LP_MPUA(LPC%ARRAY_INDEX) = B2%LP_MPUA(LPC%ARRAY_INDEX) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR_WALL)*LP%PWT*M_DROP/B1%AREA
ENDIF

! Assign particle or PARTICLE mass to the grid cell if the particle/PARTICLE not on a surface

IF (BC%IOR==0 .OR. LPC%SOLID_PARTICLE) THEN
DEN_ADD = LP%PWT*LP%MASS*LP%RVC
AREA_ADD = LP%PWT*A_DROP*LP%RVC
DROP_DEN(II,JJ,KK) = DROP_DEN(II,JJ,KK) + DEN_ADD
DROP_TMP(II,JJ,KK) = DROP_TMP(II,JJ,KK) + DEN_ADD*LP_B1%TMP_F
DROP_RAD(II,JJ,KK) = DROP_RAD(II,JJ,KK) + AREA_ADD*R_DROP
DROP_AREA(II,JJ,KK) = DROP_AREA(II,JJ,KK) + AREA_ADD
ENDIF
ENDDO PARTICLE_LOOP_2

! Compute Mass Per Unit Area (MPUA) for a liquid droplet stuck to a solid surface
! Compute cumulative quantities for PARTICLE "clouds"

IF (BC%IOR/=0 .AND. (LP%WALL_INDEX>0 .OR. LP%CFACE_INDEX>0) .AND. .NOT. SF_FIXED) THEN
IF (LP%WALL_INDEX>0) THEN
WC => WALL(LP%WALL_INDEX)
B1 => BOUNDARY_PROP1(WC%B1_INDEX)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
ELSE
CFA => CFACE(LP%CFACE_INDEX)
B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
ENDIF
IF (LPC%LIQUID_DROPLET) THEN
R_DROP = LP%RADIUS
FTPR = FOTHPI * LPC%DENSITY
M_DROP = FTPR*R_DROP**3
DO K=1,KBAR
DO J=1,JBAR
DO I=1,IBAR
DROP_RAD(I,J,K) = DROP_RAD(I,J,K)/(DROP_AREA(I,J,K)+TWO_EPSILON_EB)
DROP_TMP(I,J,K) = DROP_TMP(I,J,K)/(DROP_DEN(I,J,K) +TWO_EPSILON_EB)
AVG_DROP_RAD(I,J,K,LPC%ARRAY_INDEX ) = DROP_RAD(I,J,K)
AVG_DROP_TMP(I,J,K,LPC%ARRAY_INDEX ) = LPC%RUNNING_AVERAGE_FACTOR*AVG_DROP_TMP(I,J,K,LPC%ARRAY_INDEX ) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR)*DROP_TMP(I,J,K)
IF (LPC%Y_INDEX>0) THEN
AVG_DROP_TMP(I,J,K,LPC%ARRAY_INDEX ) = MAX(SPECIES(LPC%Y_INDEX)%TMP_MELT,AVG_DROP_TMP(I,J,K,LPC%ARRAY_INDEX ))
ELSE
M_DROP = LP%MASS
AVG_DROP_TMP(I,J,K,LPC%ARRAY_INDEX ) = MAX(TMPM,AVG_DROP_TMP(I,J,K,LPC%ARRAY_INDEX ))
ENDIF
B2%LP_MPUA(LPC%ARRAY_INDEX) = B2%LP_MPUA(LPC%ARRAY_INDEX) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR_WALL)*LP%PWT*M_DROP/B1%AREA
ENDIF

ENDDO PARTICLE_LOOP_2

! Compute cumulative quantities for PARTICLE "clouds"

DROP_RAD = DROP_RAD/(DROP_AREA+TWO_EPSILON_EB)
DROP_TMP = DROP_TMP/(DROP_DEN +TWO_EPSILON_EB)
AVG_DROP_RAD(:,:,:,LPC%ARRAY_INDEX ) = DROP_RAD
AVG_DROP_TMP(:,:,:,LPC%ARRAY_INDEX ) = LPC%RUNNING_AVERAGE_FACTOR*AVG_DROP_TMP(:,:,:,LPC%ARRAY_INDEX ) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR)*DROP_TMP
IF (LPC%Y_INDEX>0) THEN
AVG_DROP_TMP(:,:,:,LPC%ARRAY_INDEX ) = MAX(SPECIES(LPC%Y_INDEX)%TMP_MELT,AVG_DROP_TMP(:,:,:,LPC%ARRAY_INDEX ))
ELSE
AVG_DROP_TMP(:,:,:,LPC%ARRAY_INDEX ) = MAX(TMPM,AVG_DROP_TMP(:,:,:,LPC%ARRAY_INDEX ))
ENDIF
AVG_DROP_DEN(:,:,:,LPC%ARRAY_INDEX ) = LPC%RUNNING_AVERAGE_FACTOR*AVG_DROP_DEN(:,:,:,LPC%ARRAY_INDEX ) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR)*DROP_DEN
AVG_DROP_AREA(:,:,:,LPC%ARRAY_INDEX) = LPC%RUNNING_AVERAGE_FACTOR*AVG_DROP_AREA(:,:,:,LPC%ARRAY_INDEX) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR)*DROP_AREA
WHERE (AVG_DROP_DEN(:,:,:,LPC%ARRAY_INDEX )<0.0001_EB .AND. ABS(DROP_DEN)<TWO_EPSILON_EB) &
AVG_DROP_DEN(:,:,:,LPC%ARRAY_INDEX ) = 0.0_EB

! Compute mean liquid droplet temperature (LP_TEMP) on the walls

IF (LPC%LIQUID_DROPLET) THEN
AVG_DROP_DEN(I,J,K,LPC%ARRAY_INDEX ) = LPC%RUNNING_AVERAGE_FACTOR*AVG_DROP_DEN(I,J,K,LPC%ARRAY_INDEX ) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR)*DROP_DEN(I,J,K)
AVG_DROP_AREA(I,J,K,LPC%ARRAY_INDEX) = LPC%RUNNING_AVERAGE_FACTOR*AVG_DROP_AREA(I,J,K,LPC%ARRAY_INDEX) + &
(1._EB-LPC%RUNNING_AVERAGE_FACTOR)*DROP_AREA(I,J,K)
IF (AVG_DROP_DEN(I,J,K,LPC%ARRAY_INDEX )<0.0001_EB .AND. ABS(DROP_DEN(I,J,K))<TWO_EPSILON_EB) &
AVG_DROP_DEN(I,J,K,LPC%ARRAY_INDEX ) = 0.0_EB
ENDDO
ENDDO
ENDDO

! Initialize work array for storing droplet total area
! Compute mean liquid droplet temperature (LP_TEMP) on the walls

DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
WC => WALL(IW)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
B2%WORK2 = 0._EB ! drop area
ENDDO
DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
CFA => CFACE(ICF)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
B2%WORK2 = 0._EB ! drop area
ENDDO
IF (LPC%LIQUID_DROPLET) THEN

PARTICLE_LOOP_3: DO IP=1,NLP
! Initialize work array for storing droplet total area

LP => LAGRANGIAN_PARTICLE(IP)
IF (LP%CLASS_INDEX /= N_LPC) CYCLE PARTICLE_LOOP_3
BC => BOUNDARY_COORD(LP%BC_INDEX)
DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
WC => WALL(IW)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
B2%WORK2 = 0._EB ! drop area
ENDDO
DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
CFA => CFACE(ICF)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
B2%WORK2 = 0._EB ! drop area
ENDDO

IF (BC%IOR/=0 .AND. (LP%WALL_INDEX>0 .OR. LP%CFACE_INDEX>0) .AND. .NOT. SF_FIXED) THEN
IF (LP%WALL_INDEX>0) THEN
WC => WALL(LP%WALL_INDEX)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
ELSE
CFA => CFACE(LP%CFACE_INDEX)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
ENDIF
R_DROP = LP%RADIUS
A_DROP = LP%PWT*PI*R_DROP**2
B2%WORK2 = B2%WORK2 + A_DROP
ENDIF
PARTICLE_LOOP_3: DO IP=1,NLP

ENDDO PARTICLE_LOOP_3
LP => LAGRANGIAN_PARTICLE(IP)
IF (LP%CLASS_INDEX /= N_LPC) CYCLE PARTICLE_LOOP_3
BC => BOUNDARY_COORD(LP%BC_INDEX)

DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
WC => WALL(IW)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
IF (B2%WORK2 > 0._EB) THEN
B2%LP_TEMP(LPC%ARRAY_INDEX) = B2%LP_TEMP(LPC%ARRAY_INDEX)/B2%WORK2
B2%LP_TEMP(LPC%ARRAY_INDEX) = MAX(SPECIES(LPC%Y_INDEX)%TMP_MELT,B2%LP_TEMP(LPC%ARRAY_INDEX))
ELSE
B2%LP_TEMP(LPC%ARRAY_INDEX) = TMPA
ENDIF
ENDDO
DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
CFA => CFACE(ICF)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
IF (B2%WORK2 > 0._EB) THEN
B2%LP_TEMP(LPC%ARRAY_INDEX) = B2%LP_TEMP(LPC%ARRAY_INDEX)/B2%WORK2
B2%LP_TEMP(LPC%ARRAY_INDEX) = MAX(SPECIES(LPC%Y_INDEX)%TMP_MELT,B2%LP_TEMP(LPC%ARRAY_INDEX))
IF (BC%IOR/=0 .AND. (LP%WALL_INDEX>0 .OR. LP%CFACE_INDEX>0)) THEN
IF (LP%WALL_INDEX>0) THEN
WC => WALL(LP%WALL_INDEX)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
ELSE
B2%LP_TEMP(LPC%ARRAY_INDEX) = TMPA
CFA => CFACE(LP%CFACE_INDEX)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
ENDIF
ENDDO
R_DROP = LP%RADIUS
A_DROP = LP%PWT*PI*R_DROP**2
B2%WORK2 = B2%WORK2 + A_DROP
ENDIF

ENDIF
ENDDO PARTICLE_LOOP_3

DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
WC => WALL(IW)
B2 => BOUNDARY_PROP2(WC%B2_INDEX)
IF (B2%WORK2 > 0._EB) THEN
B2%LP_TEMP(LPC%ARRAY_INDEX) = B2%LP_TEMP(LPC%ARRAY_INDEX)/B2%WORK2
B2%LP_TEMP(LPC%ARRAY_INDEX) = MAX(SPECIES(LPC%Y_INDEX)%TMP_MELT,B2%LP_TEMP(LPC%ARRAY_INDEX))
ELSE
B2%LP_TEMP(LPC%ARRAY_INDEX) = TMPA
ENDIF
ENDDO
DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
CFA => CFACE(ICF)
B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
IF (B2%WORK2 > 0._EB) THEN
B2%LP_TEMP(LPC%ARRAY_INDEX) = B2%LP_TEMP(LPC%ARRAY_INDEX)/B2%WORK2
B2%LP_TEMP(LPC%ARRAY_INDEX) = MAX(SPECIES(LPC%Y_INDEX)%TMP_MELT,B2%LP_TEMP(LPC%ARRAY_INDEX))
ELSE
B2%LP_TEMP(LPC%ARRAY_INDEX) = TMPA
ENDIF
ENDDO

ENDDO PART_CLASS_SUM_LOOP
ENDIF

ENDIF SUM_PART_QUANTITIES
ENDDO PART_CLASS_SUM_LOOP

END SUBROUTINE PARTICLE_RUNNING_AVERAGES

Expand Down

0 comments on commit d8822ca

Please sign in to comment.