Merge pull request #14288 from rmcdermo/master

FDS Source: replace IMPINGIN JET model with user-specified heat trans…

Author: rmcdermo

Source/func.f90

@@ -4280,7 +4280,6 @@ SUBROUTINE PACK_BOUNDARY_PROP1(NM,IC,RC,LC,OS,B1_INDEX,UNPACK_IT,COUNT_ONLY,SURF
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_NORMAL_S,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_NORMAL_0,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_TANG,UNPACK_IT)
-RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_IMPACT,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_F,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_G,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RDN,UNPACK_IT)

Source/read.f90

@@ -7712,7 +7712,8 @@ SUBROUTINE READ_SURF(QUICK_READ)
             STRETCH_FACTOR(MAX_LAYERS),CONVECTION_LENGTH_SCALE,&
             MATL_MASS_FRACTION(MAX_LAYERS,MAX_MATERIALS),CELL_SIZE(MAX_LAYERS),CELL_SIZE_FACTOR(MAX_LAYERS),&
             EXTINCTION_TEMPERATURE,IGNITION_TEMPERATURE,HEAT_OF_VAPORIZATION,NET_HEAT_FLUX,LAYER_DIVIDE,&
-            ROUGHNESS,RADIUS,INNER_RADIUS,LENGTH,WIDTH,DT_INSERT,HEAT_TRANSFER_COEFFICIENT,HEAT_TRANSFER_COEFFICIENT_BACK,&
+            ROUGHNESS,RADIUS,INNER_RADIUS,LENGTH,WIDTH,DT_INSERT,&
+            HEAT_TRANSFER_COEFFICIENT,HEAT_TRANSFER_COEFFICIENT_BACK,HEAT_TRANSFER_COEFFICIENT_SIGMA,&
             TAU_PART,EMISSIVITY,EMISSIVITY_BACK,SPREAD_RATE,XYZ(3),MINIMUM_LAYER_THICKNESS(MAX_LAYERS),&
             VEL_GRAD,MASS_FRACTION(MAX_SPECIES),MASS_TRANSFER_COEFFICIENT,NUSSELT_C0,NUSSELT_C1,NUSSELT_C2,NUSSELT_M,&
             PARTICLE_SURFACE_DENSITY,&
@@ -7740,7 +7741,8 @@ SUBROUTINE READ_SURF(QUICK_READ)
                 DSC_CONVERSION_FACTOR,DT_INSERT,E_COEFFICIENT,&
                 EMBER_GENERATION_HEIGHT,EMBER_IGNITION_POWER_MEAN,EMBER_IGNITION_POWER_SIGMA,EMBER_TRACKING_RATIO,EMBER_YIELD,&
                 EMISSIVITY,EMISSIVITY_BACK,EXTERNAL_FLUX,EXTINCTION_TEMPERATURE,&
-                FREE_SLIP,INERT_Q_REF,FYI,GEOMETRY,HEAT_OF_VAPORIZATION,HEAT_TRANSFER_COEFFICIENT,HEAT_TRANSFER_COEFFICIENT_BACK,&
+                FREE_SLIP,INERT_Q_REF,FYI,GEOMETRY,HEAT_OF_VAPORIZATION,&
+                HEAT_TRANSFER_COEFFICIENT,HEAT_TRANSFER_COEFFICIENT_BACK,HEAT_TRANSFER_COEFFICIENT_SIGMA,&
                 HEAT_TRANSFER_MODEL,HORIZONTAL,HRRPUA,VARIABLE_THICKNESS,HT3D,ID,IGNITION_TEMPERATURE,&
                 INIT_IDS,INIT_PER_AREA,&
                 INNER_RADIUS,INTERNAL_HEAT_SOURCE,LAYER_DIVIDE,&
@@ -8322,6 +8324,13 @@ SUBROUTINE READ_SURF(QUICK_READ)
          SF%HEAT_TRANSFER_MODEL = UGENT_HTC_MODEL
    END SELECT
 
+   IF (SF%HEAT_TRANSFER_MODEL==IMPINGING_JET_HTC_MODEL .AND. &
+      (HEAT_TRANSFER_COEFFICIENT<0._EB .OR. HEAT_TRANSFER_COEFFICIENT_SIGMA<0._EB .OR. ANY(XYZ < -1.E5_EB)) ) THEN
+      WRITE(MESSAGE,'(A,A,A)') 'ERROR(XXX): SURF ',TRIM(SF%ID), &
+         ' IMPINGING JET model requires XYZ(1:3), HEAT_TRANSFER_COEFFICIENT, and HEAT_TRANSFER_COEFFICIENT_SIGMA.'
+      CALL SHUTDOWN(MESSAGE) ; RETURN
+   ENDIF
+
    SF%HRRPUA               = 1000._EB*HRRPUA
    SF%MLRPUA               = MLRPUA
    IF ((SF%HRRPUA > 0._EB .OR. SF%MLRPUA > 0) .AND. N_REACTIONS == 0) THEN
@@ -8403,6 +8412,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
    IF (HEAT_TRANSFER_COEFFICIENT_BACK < 0._EB) HEAT_TRANSFER_COEFFICIENT_BACK=HEAT_TRANSFER_COEFFICIENT
    SF%H_FIXED              = HEAT_TRANSFER_COEFFICIENT
    SF%H_FIXED_B            = HEAT_TRANSFER_COEFFICIENT_BACK
+   SF%HTC_SIGMA            = HEAT_TRANSFER_COEFFICIENT_SIGMA
    IF (RAMP_HEAT_TRANSFER_COEFFICIENT/='null') &
       CALL GET_RAMP_INDEX(RAMP_HEAT_TRANSFER_COEFFICIENT,'TIME',SF%RAMP_H_FIXED_INDEX)
    IF (RAMP_HEAT_TRANSFER_COEFFICIENT_BACK/='null') &
@@ -8758,6 +8768,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
       ENDIF
       CALL GET_RAMP_INDEX(RAMP_T_I,'T_I PROFILE',SF%RAMP_T_I_INDEX)
    ENDIF
+
    ! Boundary layer profile
 
    IF (SF%PROFILE==BOUNDARY_LAYER_PROFILE) THEN
@@ -8856,6 +8867,7 @@ SUBROUTINE SET_SURF_DEFAULTS
 HEAT_TRANSFER_MODEL     = 'null'
 HEAT_TRANSFER_COEFFICIENT = -1._EB
 HEAT_TRANSFER_COEFFICIENT_BACK = -1._EB
+HEAT_TRANSFER_COEFFICIENT_SIGMA = -1._EB
 RAMP_HEAT_TRANSFER_COEFFICIENT = 'null'
 RAMP_HEAT_TRANSFER_COEFFICIENT = 'null'
 MASS_TRANSFER_COEFFICIENT = -1._EB

Source/turb.f90

@@ -1372,7 +1372,7 @@ SUBROUTINE FORCED_CONVECTION_MODEL(NUSSELT,RE,PR_ONTH_IN,SURF_GEOMETRY_INDEX,NUS
 END SUBROUTINE FORCED_CONVECTION_MODEL
 
 
-SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_G,MU_G,VEL_G,UIMP)
+SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_G,MU_G,VEL_G)
 
 !!!!! EXPERIMENTAL !!!!!
 
@@ -1382,12 +1382,12 @@ SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_
 ! J.P. Holman, Heat Transfer, 7th Ed., McGraw-Hill, 1990, p. 346.
 
 REAL(EB), INTENT(OUT) :: H,Z_STAR
-REAL(EB), INTENT(IN) :: DZ,TMP_W,TMP_G,K_G,RHO_G,CP_G,MU_G,VEL_G,UIMP
-REAL(EB) :: NUSSELT,Q,ZC,NU_G,ALPHA_G,PR_G,D_STAR_FORCED,D_STAR_NATURAL,D_STAR_IMPACT,THETA_NATURAL,THETA_FORCED,THETA_IMPACT,&
-            Q_OLD,ERROR,DTMP,EXPON_D_IMPACT,EXPON_L_IMPACT,EXPON_T_IMPACT !,DTDN(3)
+REAL(EB), INTENT(IN) :: DZ,TMP_W,TMP_G,K_G,RHO_G,CP_G,MU_G,VEL_G
+REAL(EB) :: NUSSELT,Q,ZC,NU_G,ALPHA_G,PR_G,D_STAR_FORCED,D_STAR_NATURAL,THETA_NATURAL,THETA_FORCED,&
+            Q_OLD,ERROR,DTMP !,DTDN(3)
 INTEGER, INTENT(OUT) :: REGIME
 INTEGER :: ITER
-INTEGER, PARAMETER :: MAX_ITER=10,NATURAL=1,FORCED=2,IMPACT=3
+INTEGER, PARAMETER :: MAX_ITER=10,NATURAL=1,FORCED=2
 REAL(EB), PARAMETER :: EIGHT_NINETHS = 8._EB/9._EB
 
 REAL(EB), PARAMETER :: Z_L_NATURAL = 3.2_EB, Z_T_NATURAL=17._EB
@@ -1396,11 +1396,6 @@ SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_
 REAL(EB), PARAMETER :: Z_L_FORCED = 8._EB, Z_T_FORCED = 80._EB
 REAL(EB), PARAMETER :: C_L_FORCED = Z_L_FORCED**(-2._EB/3._EB), C_T_FORCED = C_L_FORCED * Z_T_FORCED**(2._EB/3._EB-8._EB/9._EB)
 
-REAL(EB), PARAMETER :: GAMMA_IMPACT = 1._EB
-REAL(EB), PARAMETER :: M_L_IMPACT = 0.87_EB, M_T_IMPACT = 0.87_EB
-REAL(EB), PARAMETER :: Z_L_IMPACT = 12._EB, Z_T_IMPACT = 1.E20_EB
-REAL(EB) :: C_L_IMPACT, C_T_IMPACT
-
 IF (ABS(TMP_W-TMP_G)<TWO_EPSILON_EB) THEN
    H = 0._EB
    Z_STAR = 0._EB
@@ -1415,14 +1410,6 @@ SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_
 
 THETA_NATURAL = 0.5_EB*(TMP_W+TMP_G)*K_G*ALPHA_G*NU_G/(GRAV+TWO_EPSILON_EB)
 THETA_FORCED  =     ABS(TMP_W-TMP_G)*K_G*ALPHA_G/(VEL_G+TWO_EPSILON_EB)
-THETA_IMPACT  =     ABS(TMP_W-TMP_G)*K_G*(ALPHA_G/(UIMP +TWO_EPSILON_EB))**GAMMA_IMPACT
-
-EXPON_D_IMPACT = 1._EB/(1._EB+GAMMA_IMPACT)
-EXPON_L_IMPACT = M_L_IMPACT/(1._EB+M_L_IMPACT) * (1._EB+GAMMA_IMPACT)
-EXPON_T_IMPACT = M_T_IMPACT/(1._EB+M_T_IMPACT) * (1._EB+GAMMA_IMPACT)
-
-C_L_IMPACT = Z_L_IMPACT**(-EXPON_L_IMPACT)
-C_T_IMPACT = C_L_IMPACT * Z_T_IMPACT**(EXPON_L_IMPACT-EXPON_T_IMPACT)
 
 ! ! needed if NATURAL
 ! DTDN = (TMP_G-TMP_W)/ZC * NVEC
@@ -1443,9 +1430,8 @@ SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_
    ! Step 2: compute new thermal diffusive length scale, delta*
    D_STAR_NATURAL = (THETA_NATURAL/Q)**0.25_EB
    D_STAR_FORCED  = SQRT(THETA_FORCED/Q)
-   D_STAR_IMPACT  = (THETA_IMPACT/Q)**EXPON_D_IMPACT
 
-   REGIME=MINLOC((/D_STAR_NATURAL,D_STAR_FORCED,D_STAR_IMPACT/),DIM=1)
+   REGIME=MINLOC((/D_STAR_NATURAL,D_STAR_FORCED/),DIM=1)
 
    ! Set REGIME based on minimum delta*
 
@@ -1479,20 +1465,6 @@ SUBROUTINE RAYLEIGH_HEAT_FLUX_MODEL(H,Z_STAR,REGIME,DZ,TMP_W,TMP_G,K_G,RHO_G,CP_
             NUSSELT = C_T_FORCED * Z_STAR**EIGHT_NINETHS
          ENDIF
 
-      CASE(IMPACT)
-
-         ! Step 3: compute new z* (thermal)
-         Z_STAR = ZC/D_STAR_IMPACT
-
-         ! Step 4: based on z*, choose scaling law
-         IF (Z_STAR<=Z_L_IMPACT) THEN
-            NUSSELT = 1._EB
-         ELSEIF (Z_STAR>Z_L_IMPACT .AND. Z_STAR<=Z_T_IMPACT) THEN
-            NUSSELT = C_L_IMPACT * Z_STAR**EXPON_L_IMPACT
-         ELSE
-            NUSSELT = C_T_IMPACT * Z_STAR**EXPON_T_IMPACT
-         ENDIF
-
    END SELECT REGIME_SELECT
 
    ! Step 5: update heat transfer coefficient

Source/type.f90

@@ -316,7 +316,6 @@ MODULE TYPES
    REAL(EB) :: U_NORMAL_S=0._EB      !< Estimated normal component of velocity (m/s) at next time step
    REAL(EB) :: U_NORMAL_0=0._EB      !< Initial or specified normal component of velocity (m/s) at surface
    REAL(EB) :: U_TANG=0._EB          !< Tangential velocity (m/s) near surface
-   REAL(EB) :: U_IMPACT=0._EB        !< Impact velocity from stagnation pressure
    REAL(EB) :: RHO_F                 !< Gas density at the wall (kg/m3)
    REAL(EB) :: RHO_G                 !< Gas density in near wall cell (kg/m3)
    REAL(EB) :: RDN=1._EB             !< \f$ 1/ \delta n \f$ at the surface (1/m)
@@ -879,6 +878,7 @@ MODULE TYPES
    REAL(EB) :: HM_FIXED=-1._EB                         ! Frontside mass transfer coefficient (m/s)
    REAL(EB) :: EMISSIVITY_BACK                         ! Backside emissivity
    REAL(EB) :: CONV_LENGTH                             ! Length used in heat transfer correlations
+   REAL(EB) :: HTC_SIGMA=1._EB                         ! Sigma distance of Gaussian profile in IMPINGING JET HTC model (m)
    REAL(EB) :: XYZ(3)                                  ! Starting point for a spreading fire
    REAL(EB) :: FIRE_SPREAD_RATE                        ! Defines speed (m/s) of spread from XYZ
    REAL(EB) :: INNER_RADIUS=0._EB                      ! Inner radius when SURF GEOMETRY = CYLINDRICAL
@@ -908,7 +908,7 @@ MODULE TYPES
    REAL(EB) :: M_DOT_G_PP_ADJUST_FAC=1._EB             !< For MLRPUA, scales the gas production if gas H_o_C /= solid H_o_C
    REAL(EB) :: HOC_EFF                                 !< Effective heat of combustion for S_pyro
    REAL(EB) :: Y_S_EFF                                 !< Effective soot yield for S_pyro
-   REAL(EB) :: CHI_R_EFF                              !< Effective radiative fraction for S_pyro
+   REAL(EB) :: CHI_R_EFF                               !< Effective radiative fraction for S_pyro
    REAL(EB) :: TIME_STEP_FACTOR=10._EB                 !< Maximum amount to reduce solid phase conduction time step
    REAL(EB) :: REMESH_RATIO=0.05                       !< Fraction change in wall node DX to trigger a remesh
 

Source/wall.f90

@@ -332,9 +332,6 @@ SUBROUTINE NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP,WALL_INDEX,PARTICLE_INDEX)
          VBAR = 0.5_EB*(VV(BC%IIG,BC%JJG,BC%KKG)+VV(BC%IIG,BC%JJG-1,BC%KKG)) - SF%VEL_T(2)*RAMP_FACTOR
          B1%U_TANG = SQRT(UBAR**2+VBAR**2)
    END SELECT
-   B1%U_IMPACT = SQRT(MAX(UU(BC%IIG,BC%JJG,BC%KKG),UU(BC%IIG-1,BC%JJG,BC%KKG))**2 + &
-                      MAX(VV(BC%IIG,BC%JJG,BC%KKG),VV(BC%IIG,BC%JJG-1,BC%KKG))**2 + &
-                      MAX(WW(BC%IIG,BC%JJG,BC%KKG),WW(BC%IIG,BC%JJG,BC%KKG-1))**2)
 ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
    UBAR = 0.5_EB*(UU(BC%IIG,BC%JJG,BC%KKG)+UU(BC%IIG-1,BC%JJG,BC%KKG)) - LP%U
    VBAR = 0.5_EB*(VV(BC%IIG,BC%JJG,BC%KKG)+VV(BC%IIG,BC%JJG-1,BC%KKG)) - LP%V
@@ -3515,7 +3512,7 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND
 INTEGER  :: SURF_GEOMETRY,ITMP,I,HTR
 REAL(EB) :: RE,H_NATURAL,H_FORCED,FRICTION_VELOCITY=0._EB,YPLUS=0._EB,ZSTAR,DN,TMP_FILM,MU_G,K_G,CP_G,&
             TMP_G,R_DROP,RHO_G,ZZ_G(1:N_TRACKED_SPECIES),CONV_LENGTH,GR,RA,NUSSELT_FORCED,NUSSELT_FREE,NUSSELT_IMPINGE,&
-            FILM_FAC,PHI,C0_IMP,C1_IMP,C2_IMP,M_IMP
+            FILM_FAC,PHI,XX
 INTEGER, PARAMETER :: NATURAL=1,FORCED=2,IMPACT=3,RESOLVED=4
 TYPE(MESH_TYPE), POINTER :: MX
 TYPE(SURFACE_TYPE), INTENT(IN), POINTER :: SFX
@@ -3526,7 +3523,7 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND
 TYPE(BOUNDARY_PROP2_TYPE), POINTER :: P2X
 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BCX
 
-MX  => MESHES(NMX)
+MX => MESHES(NMX)
 CONV_LENGTH = SFX%CONV_LENGTH
 
 ! Determine if this is a particle or wall cell
@@ -3583,7 +3580,7 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND
 
 ! If the user wants a specified HTC, set it and return
 
-H_FIXED_IF: IF (H_FIXED >= 0._EB) THEN
+H_FIXED_IF: IF (H_FIXED >= 0._EB .AND. SFX%HEAT_TRANSFER_MODEL/=IMPINGING_JET_HTC_MODEL) THEN
 
    HEAT_TRANSFER_COEFFICIENT = H_FIXED
 
@@ -3622,14 +3619,8 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND
             CALL NATURAL_CONVECTION_MODEL(NUSSELT_FREE,RA,SFX,BCX%IOR,DELTA_N_TMP)
             NUSSELT_IMPINGE = 0._EB
             IF (SFX%HEAT_TRANSFER_MODEL==IMPINGING_JET_HTC_MODEL) THEN
-               ! Huang, G.C. : Investigations of Heat-Transfer Coefficients for Air Flow Through Round Jets Impinging
-               ! Normal to a Heat-Transfer Surface. J. Heat Transfer, vol. 85, no. 3, Aug. 1963, pp. 237-245.
-               RE = RHO_G*P1X%U_IMPACT*CONV_LENGTH/MU_G
-               C0_IMP = 0.000_EB; IF (SFX%NUSSELT_C0>0._EB) THEN; C0_IMP = SFX%NUSSELT_C0; ENDIF
-               C1_IMP = 0.055_EB; IF (SFX%NUSSELT_C1>0._EB) THEN; C1_IMP = SFX%NUSSELT_C1; ENDIF
-               C2_IMP = 0.000_EB; IF (SFX%NUSSELT_C2>0._EB) THEN; C2_IMP = SFX%NUSSELT_C2; ENDIF
-               M_IMP  = 0.800_EB; IF (SFX%NUSSELT_M >0._EB) THEN; M_IMP  = SFX%NUSSELT_M ; ENDIF
-               NUSSELT_IMPINGE = C0_IMP + (C1_IMP*RE**M_IMP - C2_IMP) * PR_ONTH
+               XX = SQRT( (SFX%XYZ(1)-BCX%X)**2 + (SFX%XYZ(2)-BCX%Y)**2 + (SFX%XYZ(3)-BCX%Z)**2 )
+               NUSSELT_IMPINGE = HEAT_TRANSFER_COEFFICIENT*EXP(-0.5_EB * (XX / SFX%HTC_SIGMA)**2) * CONV_LENGTH/K_G
             ENDIF
             IF (PRESENT(PARTICLE_INDEX_IN)) THEN
                HEAT_TRANSFER_COEFFICIENT = MAX(NUSSELT_FORCED,NUSSELT_FREE)*K_G/CONV_LENGTH
@@ -3657,7 +3648,7 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND
          CASE(RAYLEIGH_HTC_MODEL)
             ! not appropriate for a particle, used with SURF and CFACE only
             CALL GET_SPECIFIC_HEAT(ZZ_G,CP_G,TMP_FILM)
-            CALL RAYLEIGH_HEAT_FLUX_MODEL(H_NATURAL,ZSTAR,HTR,DN,P1X%TMP_F,TMP_G,K_G,RHO_G,CP_G,MU_G,P1X%U_TANG,P1X%U_IMPACT)
+            CALL RAYLEIGH_HEAT_FLUX_MODEL(H_NATURAL,ZSTAR,HTR,DN,P1X%TMP_F,TMP_G,K_G,RHO_G,CP_G,MU_G,P1X%U_TANG)
             P2X%Z_STAR = ZSTAR
             P2X%HEAT_TRANSFER_REGIME = HTR
             HEAT_TRANSFER_COEFFICIENT = H_NATURAL