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

drjfloyd · drjfloyd · commit ba05789be637 · 2025-09-04T17:59:03.000-04:00
diff --git a/Source/devc.f90 b/Source/devc.f90
@@ -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,&
@@ -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
diff --git a/Source/dump.f90 b/Source/dump.f90
@@ -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
@@ -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
@@ -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
diff --git a/Source/read.f90 b/Source/read.f90
@@ -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,&
@@ -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
 
@@ -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
@@ -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
@@ -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)
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
