Merge pull request #15673 from drjfloyd/master

FDS Source: Add missing density factor to Bi-Dir probe.

Author: drjfloyd

Source/dump.f90

@@ -7335,13 +7335,13 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
 INTEGER, INTENT(IN) :: II,JJ,KK,IND,IND2,NM,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX, &
                        MATL_INDEX
 INTEGER, INTENT(IN), OPTIONAL :: ICC_IN,JCC_IN
-REAL(EB) :: H_TC,TMP_TC,RE_D,NUSSELT,VEL,K_G,MU_G,COSTHETA,FAC,DP,&
+REAL(EB) :: H_TC,TMP_TC,RE_D,NUSSELT,VEL,K_G,MU_G,COSTHETA,FAC,&
             Q_SUM,TMP_G,UU,VV,WW,VEL2,Y_MF_INT,PATHLENGTH,EXT_COEF,MASS_EXT_COEF,ZZ_FUEL,ZZ_OX,&
             VELSR,WATER_VOL_FRAC,RHS,DT_C,DT_E,T_RATIO,Y_E_LAG, H_G,H_G_SUM,CPBAR,CP,ZZ_GET(1:N_TRACKED_SPECIES),RCON,&
             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,PROBE_TMP
+            GAS_PHASE_OUTPUT_CFA,CFACE_AREA,VELOCITY_COMPONENT(1:3),ATOTV(1:3),TMP_F,R_D,MW,PROBE_TMP,PROBE_RHO
 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
@@ -7893,28 +7893,33 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
       GAS_PHASE_OUTPUT_RES = (GAS_PHASE_OUTPUT_RES - K_G*(TMP(IP,JP,KP)-TMP(II,JJ,KK))*R_DN)*0.001
 
    CASE(114) ! BI-DIRECTIONAL PROBE
+      ! Fits taken from
+      ! McCaffrey and Heskestad, A Robust Bidirectional Low-Velocity Probe for Flame and Fire Application
+      ! Combustion and Flame, 26, 125 - 127, (1976).
       IF (PY%TC) THEN
          PROBE_TMP = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,110,IND2,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,&
                                            PROP_INDEX,REAC_INDEX,MATL_INDEX,ICC_IN,JCC_IN) + TMPM
       ELSE
          PROBE_TMP = TMP(II,JJ,KK)
       ENDIF
+      PROBE_RHO = MW_AIR*P_STP/(R0*PROBE_TMP)
       UU = 0.5_EB*(U(MAX(0,II-1),JJ,KK)+U(MIN(IBAR,II),JJ,KK))
       VV = 0.5_EB*(V(II,MAX(0,JJ-1),KK)+V(II,MIN(JBAR,JJ),KK))
       WW = 0.5_EB*(W(II,JJ,MAX(0,KK-1))+W(II,JJ,MIN(KBAR,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)+&
+      ! Adjust for effect of flow direction on measured pressure
+      COSTHETA = (UU*ORIENTATION_VECTOR(1,DV%ORIENTATION_INDEX)+VV*ORIENTATION_VECTOR(2,DV%ORIENTATION_INDEX)+ &
                   WW*ORIENTATION_VECTOR(3,DV%ORIENTATION_INDEX))/VEL
+      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
+      ! Adjust for effect of Re number on measured pressure
       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
-      GAS_PHASE_OUTPUT_RES = SIGN(1._EB,COSTHETA)*VEL*PY%CALIBRATION_CONSTANT*FAC
+      GAS_PHASE_OUTPUT_RES = SIGN(1._EB,COSTHETA)*VEL*PY%CALIBRATION_CONSTANT*FAC*SQRT(RHO(II,JJ,KK)/PROBE_RHO)
 
    CASE(130) ! EXTINCTION
       ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)