FDS Source: Allow for dynamic bi-dir probe calibration constant

Author: drjfloyd

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -10567,14 +10567,14 @@ \subsection{Bi-Directional Probe}
 \label{info:bidir_probe}
 \label{bi_dir}
 
-The output quantity \ct{BI-DIRECTIONAL PROBE} is the velocity of a modeled bi-directional probe. A bi-directional probe uses the following equation:
+The output quantity \ct{BI-DIRECTIONAL PROBE} is the velocity of a modeled bi-directional probe. A bi-directional probe uses the following equation~\cite{McCaffrey:1976}:
 \be
 C \sqrt{\frac{2 \Delta P}{\rho}}
 \label{BDP}
 \ee
 where $C$ is a calibration constant (default value is 0.93), $\Delta P$ is the pressure difference across the probe, and $\rho$ is the gas density at the probe. In a typical experiment, the gas density is computed assuming standard pressure (101325 Pa), the molecular weight of air (28.8 g/mol), and the temperature as measured by a thermocouple near the probe.
 
-Bi-directional probes have biases due to both the Reynolds number (based on the probe diameter) of the flow and the angle of the flow with respect to the probe axis~\cite{McCaffrey:1976}. At low Reynolds number a probe will measure a higher effective velocity. As the angle of the flow vector with the axis increases, the effective velocity at first increases up to an angle of 30$^\circ$ due to a low pressure region forming downstream of the probe, and then decreases reaching no measured flow at an angle of 90$^\circ$. This model accounts for these sensitivities and the impact of density differences from varied molecular weight at the probe. The orientation of the probe can be specified with either \ct{IOR} or \ct{ORIENTATION} on \ct{DEVC}. A probe with \ct{IOR}=-1 would have a positive velocity output when the flow is in the -x direction. Parameters for the probe can be specified with a \ct{PROP_ID} on the \ct {DEVC}. The calibration constant (default of 0.93) and the probe diameter (default of 0.0254 m) can be set respectively with \ct{CALIBRATION_CONSTANT} and \ct{PROBE_DIAMETER} on \ct{PROP}. If the probe temperature is an aspirated thermocouple or other measurement not sensitive to the radiative environment, then set \ct{TC=F} on \ct{PROP}. Thermocouple specific properties for a bi-directional probe, see Section~\ref{info:THERMOCOUPLE}, should be set with the same \ct{PROP} as for the probe.
+Bi-directional probes have biases due to both the Reynolds number (based on the probe diameter) of the flow and the angle of the flow with respect to the probe axis~\cite{McCaffrey:1976}. At low Reynolds number a probe will measure a higher effective velocity. As the angle of the flow vector with the axis increases, the effective velocity at first increases up to an angle of 30$^\circ$ due to a low pressure region forming downstream of the probe, and then decreases reaching no measured flow at an angle of 90$^\circ$. This model accounts for these sensitivities and the impact of density differences from varied molecular weight at the probe. The orientation of the probe can be specified with either \ct{IOR} or \ct{ORIENTATION} on \ct{DEVC}. A probe with \ct{IOR}=-1 would have a positive velocity output when the flow is in the -x direction. Parameters for the probe can be specified with a \ct{PROP_ID} on the \ct {DEVC}. The calibration constant (default of 0.93) and the probe diameter (default of 0.0254 m) can be set respectively with \ct{CALIBRATION_CONSTANT} and \ct{PROBE_DIAMETER} on \ct{PROP}. If the probe temperature is an aspirated thermocouple or other measurement not sensitive to the radiative environment, then set \ct{TC=F} on \ct{PROP}. Thermocouple specific properties for a bi-directional probe, see Section~\ref{info:THERMOCOUPLE}, should be set with the same \ct{PROP} as for the probe. If a dynamic calibration constant based on the Reynolds number calibration in McCaffrey~\cite{McCaffrey:1976} was used set \ct{CALIBRATION_CONSTANT=-1}.
 
 Figure~\ref{bi_dir_fig} shows the results of a bi-directional probe with varying angle to a 1~m/s flow and varying flow speed.
 \begin{figure}[ht]

Source/dump.f90

@@ -7341,7 +7341,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,PROBE_TMP,PROBE_RHO
+            GAS_PHASE_OUTPUT_CFA,CFACE_AREA,VELOCITY_COMPONENT(1:3),ATOTV(1:3),TMP_F,R_D,MW,PROBE_TMP,PROBE_RHO,PROBE_DELTA_P
 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
@@ -7918,8 +7918,19 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
       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 = 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*SQRT(RHO(II,JJ,KK)/PROBE_RHO)
+            0.0000000000002555_EB*RE_D**4-2.484E-17_EB*RE_D**5
+      IF (PY%CALIBRATION_CONSTANT > 0._EB) THEN
+         GAS_PHASE_OUTPUT_RES = SIGN(1._EB,COSTHETA)*VEL*PY%CALIBRATION_CONSTANT*FAC*SQRT(RHO(II,JJ,KK)/PROBE_RHO)
+      ELSE
+         PROBE_DELTA_P = (VEL*FAC)**2*RHO(II,JJ,KK)*0.5_EB
+         ! LJ AIR viscosity fit 100 K to 5000 K
+         MU_G = = 1.5205E-22_EB*PROBE_TMP**5 - 2.1417E-18_EB*PROBE_TMP**4 + 1.1402E-14_EB*PROBE_TMP**3 - &
+                  2.9846E-11_EB*PROBE_TMP**2 + 5.9898E-8_EB*PROBE_TMP + 0.000002352_EB
+         RE_D = MIN(3800._EB,MAX(40._EB,PROBE_RHO*VEL*PY%PROBE_DIAMETER/MU_G))
+         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)*1._EB/FAC*SQRT(2*PROBE_DELTA_P/PROBE_RHO)
+      ENDIF
 
    CASE(130) ! EXTINCTION
       ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)