Skip to content

FDS Source: BI-DIRECTIONAL PROBE edits and Type-K props #15064

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

Merged
drjfloyd merged 1 commit into from
Sep 4, 2025
Merged

FDS Source: BI-DIRECTIONAL PROBE edits and Type-K props #15064

Show file tree
Hide file tree
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
7 changes: 4 additions & 3 deletions Source/devc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ MODULE DEVICE_VARIABLES
!> \brief Derived type storing inputs for the PROP namelist group

TYPE PROPERTY_TYPE
REAL(EB) :: DENSITY,DIAMETER,EMISSIVITY,HEAT_TRANSFER_COEFFICIENT,SPECIFIC_HEAT,RTI,TIME_CONSTANT, &
REAL(EB) :: DENSITY,DIAMETER,EMISSIVITY,HEAT_TRANSFER_COEFFICIENT,RTI,TIME_CONSTANT, &
ACTIVATION_TEMPERATURE,ACTIVATION_OBSCURATION, CALIBRATION_CONSTANT,&
ALPHA_E,ALPHA_C,BETA_E,BETA_C,CHARACTERISTIC_VELOCITY,PARTICLE_VELOCITY,MASS_FLOW_RATE,FLOW_RATE,FLOW_TAU, &
GAUGE_EMISSIVITY,GAUGE_TEMPERATURE,INITIAL_TEMPERATURE,K_FACTOR,C_FACTOR,OPERATING_PRESSURE,OFFSET,&
Expand All @@ -19,8 +19,9 @@ MODULE DEVICE_VARIABLES
LOGICAL :: PDPA_INTEGRATE=.TRUE.,PDPA_NORMALIZE=.TRUE.,HISTOGRAM_NORMALIZE=.TRUE.,HISTOGRAM=.FALSE., &
HISTOGRAM_CUMULATIVE=.FALSE.,SPARK=.FALSE.,TC=.TRUE.,IGNITION_ZONE=.FALSE.
REAL(EB) :: PDPA_START=0._EB,PDPA_END=1.E6_EB,PDPA_RADIUS=0.1_EB
REAL(EB), ALLOCATABLE, DIMENSION(:) :: TABLE_ROW, V_FACTOR
INTEGER :: PART_INDEX=-1,FLOW_RAMP_INDEX,SPRAY_PATTERN_INDEX,Z_INDEX=-999,Y_INDEX=-999,PRESSURE_RAMP_INDEX
REAL(EB), ALLOCATABLE, DIMENSION(:) :: TABLE_ROW, V_FACTOR,SPECIFIC_HEAT
INTEGER :: PART_INDEX=-1,FLOW_RAMP_INDEX,SPRAY_PATTERN_INDEX,Z_INDEX=-999,Y_INDEX=-999,PRESSURE_RAMP_INDEX,&
SPECIFIC_HEAT_RAMP_INDEX=-1
CHARACTER(LABEL_LENGTH) :: SMOKEVIEW_ID(SMOKEVIEW_OBJECTS_DIMENSION),PART_ID,ID,QUANTITY,TABLE_ID,SPEC_ID='null', &
SMOKEVIEW_PARAMETERS(SMOKEVIEW_OBJECTS_DIMENSION)='null'
REAL(EB), ALLOCATABLE, DIMENSION(:) :: SPRAY_LON_CDF,SPRAY_LON,SPRAY_LAT
Expand Down
47 changes: 25 additions & 22 deletions Source/dump.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7334,7 +7334,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
EXPON,Y_SPECIES,MEC,Y_SPECIES2,Y_H2O,R_Y_H2O,R_DN,SGN,Y_ALL(N_SPECIES),H_S,D_Z_N(0:I_MAX_TEMP),&
DISSIPATION_RATE,S11,S22,S33,S12,S13,S23,DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ,ONTHDIV,SS,ETA,DELTA,R_DX2,&
UVW,UODX,VODY,WODZ,XHAT,ZHAT,BBF,GAMMA_LOC,VC,VOL,PHI,GAS_PHASE_OUTPUT_CC,&
GAS_PHASE_OUTPUT_CFA,CFACE_AREA,VELOCITY_COMPONENT(1:3),ATOTV(1:3),TMP_F,R_D,MW,RHO_AIR,PROBE_TMP
GAS_PHASE_OUTPUT_CFA,CFACE_AREA,VELOCITY_COMPONENT(1:3),ATOTV(1:3),TMP_F,R_D,MW,PROBE_TMP
INTEGER :: N,I,J,K,NN,IL,III,JJJ,KKK,IP,JP,KP,FED_ACTIVITY,IP1,JP1,KP1,IM1,JM1,KM1,IIM1,JJM1,KKM1,NR,NS,RAM,&
ICC,JCC,NCELL,AXIS,ICF,NFACE,JCF,JCC_LO,JCC_HI,PDPA_FORMULA,IC
REAL(FB) :: RN
Expand Down Expand Up @@ -7843,24 +7843,27 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
GAS_PHASE_OUTPUT_RES = FED(ZZ_GET,RSUM(II,JJ,KK),FED_ACTIVITY)

CASE(110) ! THERMOCOUPLE
TMP_G = TMP(II,JJ,KK)
IF (PY%HEAT_TRANSFER_COEFFICIENT<0._EB) THEN
UU = U(II,JJ,KK)
VV = V(II,JJ,KK)
WW = W(II,JJ,KK)
VEL2 = UU**2+VV**2+WW**2
ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP(II,JJ,KK))
CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP(II,JJ,KK))
RE_D = RHO(II,JJ,KK)*SQRT(VEL2)*PY%DIAMETER/MU_G
CALL FORCED_CONVECTION_MODEL(NUSSELT,RE_D,PR_ONTH,SURF_SPHERICAL)
H_TC = NUSSELT*K_G/PY%DIAMETER
ELSE
H_TC = PY%HEAT_TRANSFER_COEFFICIENT
IF (T > T_BEGIN) THEN
TMP_G = TMP(II,JJ,KK)
IF (PY%HEAT_TRANSFER_COEFFICIENT<0._EB) THEN
UU = U(II,JJ,KK)
VV = V(II,JJ,KK)
WW = W(II,JJ,KK)
VEL2 = UU**2+VV**2+WW**2
ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP(II,JJ,KK))
CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP(II,JJ,KK))
RE_D = RHO(II,JJ,KK)*SQRT(VEL2)*PY%DIAMETER/MU_G
CALL FORCED_CONVECTION_MODEL(NUSSELT,RE_D,PR_ONTH,SURF_SPHERICAL)
H_TC = NUSSELT*K_G/PY%DIAMETER
ELSE
H_TC = PY%HEAT_TRANSFER_COEFFICIENT
ENDIF
FAC = 6._EB/(PY%DENSITY*PY%SPECIFIC_HEAT(NINT(DV%TMP_L))*PY%DIAMETER)*DT
DV%TMP_L = (DV%TMP_L + FAC*(H_TC*(TMP_G-0.5_EB*DV%TMP_L) + &
PY%EMISSIVITY*(0.25_EB*UII(II,JJ,KK)+SIGMA*DV%TMP_L**4))) / &
(1._EB + FAC*(0.5_EB*H_TC+2._EB*PY%EMISSIVITY*SIGMA*DV%TMP_L**3))
ENDIF
RHS = (6._EB/(PY%DENSITY*PY%SPECIFIC_HEAT*PY%DIAMETER))* &
( H_TC*(TMP_G-DV%TMP_L) + PY%EMISSIVITY*(0.25_EB*UII(II,JJ,KK)-SIGMA*DV%TMP_L**4) )
IF (T>T_BEGIN) DV%TMP_L = DV%TMP_L + DT*RHS
GAS_PHASE_OUTPUT_RES = DV%TMP_L - TMPM

CASE(111:113) ! ENTHALPY FLUX
Expand Down Expand Up @@ -7893,17 +7896,17 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
VV = V(II,JJ,KK)
WW = W(II,JJ,KK)
VEL2 = UU**2+VV**2+WW**2
VEL = SQRT(VEL2)
DP = 0.5_EB*VEL2*RHO(II,JJ,KK)
COSTHETA = (UU*ORIENTATION_VECTOR(1,DV%ORIENTATION_INDEX)+VV*ORIENTATION_VECTOR(2,DV%ORIENTATION_INDEX)+&
WW*ORIENTATION_VECTOR(3,DV%ORIENTATION_INDEX))/SQRT(VEL2)
FAC = -2.308_EB*ABS(COSTHETA)**3 + 2.533_EB*ABS(COSTHETA)**2 + 0.7847_EB*ABS(COSTHETA)
VEL = FAC*SQRT(VEL2)
WW*ORIENTATION_VECTOR(3,DV%ORIENTATION_INDEX))/VEL
ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP(II,JJ,KK))
RE_D = MIN(3800._EB,MAX(40._EB,RHO(II,JJ,KK)*VEL*PY%PROBE_DIAMETER/MU_G))
FAC = MAX(0._EB,-2.308_EB*ABS(COSTHETA)**3 + 2.533_EB*ABS(COSTHETA)**2 + 0.7847_EB*ABS(COSTHETA) - 0.0097_EB)
VEL = FAC*VEL
FAC = 1.533_EB-0.001366_EB*RE_D+0.000001688_EB*RE_D**2-0.0000000009706_EB*RE_D**3+&
0.0000000000002555_EB*RE_D**4-2.484E-17_EB*RE_D**5
RHO_AIR = 350.9736_EB/PROBE_TMP !350 is 0.0288 101325/ 8.314472
GAS_PHASE_OUTPUT_RES = SIGN(1._EB,COSTHETA)*VEL*PY%CALIBRATION_CONSTANT*FAC

CASE(130) ! EXTINCTION
Expand Down
52 changes: 41 additions & 11 deletions Source/read.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6540,7 +6540,7 @@ SUBROUTINE READ_PROP
CHARACTER(LABEL_LENGTH) :: SMOKEVIEW_ID(SMOKEVIEW_OBJECTS_DIMENSION),QUANTITY='null',PART_ID='null',FLOW_RAMP='null', &
SPRAY_PATTERN_TABLE='null',SPEC_ID='null',&
PRESSURE_RAMP='null',SMOKEVIEW_PARAMETERS(SMOKEVIEW_OBJECTS_DIMENSION), &
SPRAY_PATTERN_SHAPE='GAUSSIAN'
SPRAY_PATTERN_SHAPE='GAUSSIAN',SPECIFIC_HEAT_RAMP
TYPE (PROPERTY_TYPE), POINTER :: PY

NAMELIST /PROP/ ACTIVATION_OBSCURATION,ACTIVATION_TEMPERATURE,ALPHA_C,ALPHA_E,BETA_C,BETA_E,CALIBRATION_CONSTANT,&
Expand All @@ -6551,7 +6551,7 @@ SUBROUTINE READ_PROP
PARTICLES_PER_SECOND,PARTICLE_VELOCITY,PART_ID,PDPA_END,&
PDPA_INTEGRATE,PDPA_M,PDPA_N,PDPA_NORMALIZE,PDPA_RADIUS,&
PDPA_START,PRESSURE_RAMP,PROBE_DIAMETER,PX,PXX,QUANTITY,RTI,SMOKEVIEW_ID,SMOKEVIEW_PARAMETERS,SPARK,&
SPEC_ID,SPECIFIC_HEAT,SPRAY_ANGLE,&
SPEC_ID,SPECIFIC_HEAT,SPECIFIC_HEAT_RAMP,SPRAY_ANGLE,&
SPRAY_PATTERN_BETA,SPRAY_PATTERN_MU,SPRAY_PATTERN_SHAPE,SPRAY_PATTERN_TABLE,TC,TIME_CONSTANT,VELOCITY_COMPONENT,&
VIEW_ANGLE

Expand Down Expand Up @@ -6600,7 +6600,16 @@ SUBROUTINE READ_PROP
PY%DIAMETER = DIAMETER
PY%EMISSIVITY = EMISSIVITY
PY%HEAT_TRANSFER_COEFFICIENT= HEAT_TRANSFER_COEFFICIENT
PY%SPECIFIC_HEAT = SPECIFIC_HEAT*1000._EB/TIME_SHRINK_FACTOR
ALLOCATE(PY%SPECIFIC_HEAT(0:I_MAX_TEMP))
IF(SPECIFIC_HEAT > 0._EB) THEN
PY%SPECIFIC_HEAT = SPECIFIC_HEAT*1000._EB/TIME_SHRINK_FACTOR
ELSE
! Type-K CP(20 C)=0.4515, CP(1200 C)=0.6010
DO J = 0,I_MAX_TEMP
PY%SPECIFIC_HEAT(J) = 0.4515_EB+0.001_EB*(REAL(MAX(20,MIN(1200,J)),EB)-20._EB)*(0.6010_EB-0.4515_EB)
ENDDO
ENDIF
IF (SPECIFIC_HEAT_RAMP /= 'null') CALL GET_RAMP_INDEX(SPECIFIC_HEAT_RAMP,'TEMPERATURE',PY%SPECIFIC_HEAT_RAMP_INDEX)
PY%C_FACTOR = C_FACTOR
PY%CHARACTERISTIC_VELOCITY = CHARACTERISTIC_VELOCITY
PY%GAUGE_EMISSIVITY = GAUGE_EMISSIVITY
Expand Down Expand Up @@ -6842,11 +6851,12 @@ SUBROUTINE SET_PROP_DEFAULTS
BETA_C = -1.0_EB
BETA_E = -1.0_EB
CALIBRATION_CONSTANT = 0.93_EB
DENSITY = 8908._EB ! kg/m3 (Nickel)
DENSITY = 8700._EB ! kg/m3 (Type-K)
DIAMETER = 0.001 ! m
EMISSIVITY = 0.85_EB
HEAT_TRANSFER_COEFFICIENT= -1._EB ! W/m2/K
SPECIFIC_HEAT = 0.44_EB ! kJ/kg/K (Nickel)
SPECIFIC_HEAT = -1._EB ! kJ/kg/K
SPECIFIC_HEAT_RAMP = 'null'
C_FACTOR = 0.0_EB
CHARACTERISTIC_VELOCITY = 1.0_EB ! m/s
PARTICLE_VELOCITY = 0._EB ! m/s
Expand Down Expand Up @@ -6913,6 +6923,7 @@ END SUBROUTINE READ_PROP
SUBROUTINE PROC_PROP

USE DEVICE_VARIABLES
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
REAL(EB) :: TOTAL_FLOWRATE, SUBTOTAL_FLOWRATE
INTEGER :: N,NN,N_V_FACTORS,ILPC
LOGICAL :: TABLE_NORMED(1:N_TABLE)
Expand Down Expand Up @@ -7011,7 +7022,21 @@ SUBROUTINE PROC_PROP
PY%V_FACTOR = PY%PARTICLE_VELOCITY/SQRT(PY%OPERATING_PRESSURE)
ENDIF
ENDIF

! Check units of specific heat

IF (PY%SPECIFIC_HEAT_RAMP_INDEX > 0) THEN
IF (.NOT.RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%DEP_VAR_UNITS_CONVERTED) THEN
RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%INTERPOLATED_DATA(:) = &
RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%INTERPOLATED_DATA(:)*1000._EB/TIME_SHRINK_FACTOR
RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%DEP_VAR_UNITS_CONVERTED = .TRUE.
ENDIF
IF (RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%DEPENDENT_DATA(1) > 10._EB .AND. MY_RANK==0) &
WRITE(LU_ERR,'(A,A)') 'WARNING: SPECIFIC_HEAT units are kJ/kg/K check PROP ',TRIM(PY%ID)
DO I=0,I_MAX_TEMP
PY%SPECIFIC_HEAT(I)=EVALUATE_RAMP(REAL(I,EB),PY%SPECIFIC_HEAT_RAMP_INDEX)
ENDDO
ENDIF
ENDDO PROP_LOOP

END SUBROUTINE PROC_PROP
Expand Down Expand Up @@ -14161,8 +14186,7 @@ SUBROUTINE READ_DEVC
DV%QUANTITY(1)=='GAUGE HEAT FLUX GAS' .OR. &
DV%QUANTITY(1)=='RADIANCE' .OR. &
DV%QUANTITY(1)=='ADIABATIC SURFACE TEMPERATURE GAS' .OR. &
DV%QUANTITY(1)=='RADIOMETER GAS' .OR. &
DV%QUANTITY(1)=='BI-DIRECTIONAL PROBE') THEN
DV%QUANTITY(1)=='RADIOMETER GAS') THEN
IF (DV%ORIENTATION_INDEX==0) THEN
WRITE(MESSAGE,'(3A)') 'ERROR(887): DEVC ',TRIM(ID),' must have an ORIENTATION.'
CALL SHUTDOWN(MESSAGE) ; RETURN
Expand All @@ -14175,6 +14199,12 @@ SUBROUTINE READ_DEVC
INIT_RESERVED(N_INIT_RESERVED)%DEVC_INDEX = N_DEVC
INIT_RESERVED(N_INIT_RESERVED)%N_PARTICLES = POINTS
ENDIF
ENDIF

! Special case for QUANTITY='BI-DIRECTIONAL PROBE'
IF (DV%QUANTITY(1)=='BI-DIRECTIONAL PROBE' .AND. DV%ORIENTATION_INDEX==0) THEN
WRITE(MESSAGE,'(3A)') 'ERROR(887): DEVC ',TRIM(ID),' must have an ORIENTATION.'
CALL SHUTDOWN(MESSAGE) ; RETURN
ENDIF

! Miscellaneous actions taken based on specific device attributes
Expand Down Expand Up @@ -15161,7 +15191,7 @@ SUBROUTINE PROC_DEVC
PY => PROPERTY(DV%PROP_INDEX)
IF (PY%TIME_CONSTANT>0._EB) THEN ! Convert a specified TIME_CONSTANT into an equivalent bead DIAMETER
IF (PY%HEAT_TRANSFER_COEFFICIENT>0._EB) THEN ! Calculate effective diameter directly
PY%DIAMETER = 6._EB*PY%HEAT_TRANSFER_COEFFICIENT*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT)
PY%DIAMETER = 6._EB*PY%HEAT_TRANSFER_COEFFICIENT*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT(293))
ELSE ! Calculate effective diameter implicitly
ZZ_G = 0._EB
ZZ_G(1) = 1._EB
Expand All @@ -15171,13 +15201,13 @@ SUBROUTINE PROC_DEVC
TOL = 1._EB
RE = RHOA*UU*PY%DIAMETER/MU_G ! Make initial estimate of Re and Nu based on default bead diameter
NU = 2._EB + 0.6_EB*SQRT(RE)*PR_AIR**ONTH
PY%DIAMETER = SQRT(6._EB*K_G*NU*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT))
PY%DIAMETER = SQRT(6._EB*K_G*NU*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT(293)))
DO WHILE(TOL>1.E-5_EB)
RE = RHOA*UU*PY%DIAMETER/MU_G
NU = 2._EB + 0.6_EB*SQRT(RE)*PR_AIR**ONTH
F = 6._EB*K_G*NU*PY%TIME_CONSTANT - PY%DENSITY*PY%SPECIFIC_HEAT*PY%DIAMETER**2
F = 6._EB*K_G*NU*PY%TIME_CONSTANT - PY%DENSITY*PY%SPECIFIC_HEAT(293)*PY%DIAMETER**2
DFDD = 1.8_EB*K_G*PY%TIME_CONSTANT*PR_AIR**ONTH*RE**(-0.5)*RHOA*UU/MU_G - &
2._EB*PY%DENSITY*PY%SPECIFIC_HEAT*PY%DIAMETER
2._EB*PY%DENSITY*PY%SPECIFIC_HEAT(293)*PY%DIAMETER
PY%DIAMETER = PY%DIAMETER - F/DFDD
TOL = ABS(F/DFDD)
ENDDO
Expand Down
Loading