Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FDS Source: Remove SF_FIXED from PARTICLE_RUNNING_AVERAGES and do som… #12225

Merged
merged 1 commit into from
Nov 12, 2023
Merged

FDS Source: Remove SF_FIXED from PARTICLE_RUNNING_AVERAGES and do som… #12225

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
285 changes: 139 additions & 146 deletions Source/part.f90
View file
Open in desktop
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