diff --git a/Build/impi_intel_linux_openmp_db/README.md b/Build/impi_intel_linux_openmp_db/README.md new file mode 100644 index 00000000000..ea7f662b0eb --- /dev/null +++ b/Build/impi_intel_linux_openmp_db/README.md @@ -0,0 +1,11 @@ +## Detecting OpenMP Thread Race Conditions + +To check for OpenMP race conditions, add `-fsanitize=thread` to the arguments (`FFLAGS`) for `impi_intel_linux_openmp_db` and add the lines +``` +export OMP_TOOL_LIBRARIES='libarcher.so' +export TSAN_OPTIONS='ignore_noninstrumented_modules=1' +``` +to the bash script containing the PBS/Torque or Slurm directives. + +Then run a short, small job and race errors should be written to the `.err` file. + diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 098f72c9c3c..ca99922f208 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -2299,7 +2299,7 @@ \subsubsection{Output for Convective Heat Transfer Regime} \subsubsection{Specified Convective Heat Transfer Coefficient} -To specify the convective heat transfer coefficient, set it to a constant using \ct{HEAT_TRANSFER_COEFFICIENT} on the \ct{SURF} line in units of \unit{W/(m^2.K)} with optional time dependent ramp using \ct{RAMP_HEAT_TRANSFER_COEFFICIENT}. If the back side of the solid obstruction faces the exterior of the computational domain and the solid conducts heat, the heat transfer coefficient of the back side may be specified using \ct{HEAT_TRANSFER_COEFFICIENT_BACK} with optional time dependent ramp ramp using \ct{RAMP_HEAT_TRANSFER_COEFFICIENT_BACK}. This back side condition is appropriate for a \ct{SURF} line with \ct{BACKING='VOID'} or \ct{BACKING='EXPOSED'}. +To specify the convective heat transfer coefficient, set it to a constant using \ct{HEAT_TRANSFER_COEFFICIENT} on the \ct{SURF} line in units of \unit{W/(m^2.K)} with optional time dependent ramp using \ct{RAMP_HEAT_TRANSFER_COEFFICIENT}. If the back side of the solid obstruction faces the exterior of the computational domain and the solid conducts heat, the heat transfer coefficient of the back side may be specified using \ct{HEAT_TRANSFER_COEFFICIENT_BACK} with optional time dependent ramp ramp using \ct{RAMP_HEAT_TRANSFER_COEFFICIENT_BACK}. This back side condition is appropriate for a \ct{SURF} line with \ct{BACKING='VOID'} or if the solid obsttuction backs to the exterior of the computational domain. \subsubsection{Impinging Jet Heat Transfer Model} \label{info:impinging_jet} diff --git a/Source/cons.f90 b/Source/cons.f90 index 3372871422d..43c297ed42f 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -346,8 +346,8 @@ MODULE GLOBAL_CONSTANTS REAL(EB), ALLOCATABLE, DIMENSION(:) :: H_V_H2O !< Heat of vaporization for water (J/kg) REAL(EB) :: CHI_R_MIN=0._EB !< Lower bound for radiative fraction REAL(EB) :: CHI_R_MAX=1._EB !< Upper bound for radiative fraction -REAL(EB) :: SPHERE_FILM_FAC=ONTH !< Factor used in droplet evaporation algorithm for droplets -REAL(EB) :: PLATE_FILM_FAC=ONTH !< Factor used in droplet evaporation algorithm for walls +REAL(EB) :: SPHERE_FILM_FACTOR=ONTH !< Weighting factor used in droplet evaporation algorithm for droplets +REAL(EB) :: PLATE_FILM_FACTOR=ONTH !< Weighting factor used in droplet evaporation algorithm for walls REAL(EB) :: ORIGIN_LAT=-1.E6_EB !< Latitude of terrain map REAL(EB) :: ORIGIN_LON=-1.E6_EB !< Longitude of terrain map REAL(EB) :: NORTH_BEARING=0._EB !< North bearing for terrain map diff --git a/Source/part.f90 b/Source/part.f90 index 9b94c8c3645..04414d61960 100644 --- a/Source/part.f90 +++ b/Source/part.f90 @@ -3946,7 +3946,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) SOLID_OR_GAS_PHASE_2: IF (BC%IOR/=0 .AND. (LP%WALL_INDEX>0 .OR. LP%CFACE_INDEX>0) .AND. .NOT. SF_FIXED) THEN - CALL GET_FILM_PROPERTIES(1,PLATE_FILM_FAC,Y_DROP_A,Y_GAS_A,Z_INDEX_A,TMP_DROP,TMP_G,ZZ_GET, & + CALL GET_FILM_PROPERTIES(1,PLATE_FILM_FACTOR,Y_DROP_A,Y_GAS_A,Z_INDEX_A,TMP_DROP,TMP_G,ZZ_GET, & PBAR(KK,PRESSURE_ZONE(II,JJ,KK)),TMP_FILM,MU_FILM,& K_FILM,CP_FILM,D_FILM,RHO_FILM,PR_FILM,SC_FILM) @@ -4018,7 +4018,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) ELSE SOLID_OR_GAS_PHASE_2 - CALL GET_FILM_PROPERTIES(1,SPHERE_FILM_FAC,Y_DROP_A,Y_GAS_A,Z_INDEX_A,TMP_DROP,TMP_G,ZZ_GET,& + CALL GET_FILM_PROPERTIES(1,SPHERE_FILM_FACTOR,Y_DROP_A,Y_GAS_A,Z_INDEX_A,TMP_DROP,TMP_G,ZZ_GET,& PBAR(KK,PRESSURE_ZONE(II,JJ,KK)),TMP_FILM,MU_FILM,& K_FILM,CP_FILM,D_FILM,RHO_FILM,PR_FILM,SC_FILM) diff --git a/Source/read.f90 b/Source/read.f90 index e570556ebe4..079dbea38e9 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -8305,6 +8305,13 @@ SUBROUTINE READ_SURF(QUICK_READ) CALL SHUTDOWN(MESSAGE) ; RETURN END SELECT + SELECT CASE(SF%GEOMETRY) + CASE DEFAULT ; SF%FILM_FACTOR = PLATE_FILM_FACTOR + CASE(SURF_SPHERICAL) ; SF%FILM_FACTOR = SPHERE_FILM_FACTOR + CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; SF%FILM_FACTOR = PLATE_FILM_FACTOR + CASE(SURF_CARTESIAN) ; SF%FILM_FACTOR = PLATE_FILM_FACTOR + END SELECT + SF%H_V = 1000._EB*HEAT_OF_VAPORIZATION SELECT CASE(HEAT_TRANSFER_MODEL) CASE DEFAULT @@ -9517,6 +9524,7 @@ SUBROUTINE PROC_SURF_2 SF%SPECIES_BC_INDEX = HVAC_BOUNDARY ENDIF IF (N==MASSLESS_TARGET_SURF_INDEX) THEN + SF%H_FIXED = 10._EB SF%EMISSIVITY = 1._EB ENDIF diff --git a/Source/type.f90 b/Source/type.f90 index c7d928367dd..68c62a8c210 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -910,6 +910,7 @@ MODULE TYPES 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 + REAL(EB) :: FILM_FACTOR=ONTH !< Weighting factor for evaluating surface file properties REAL(EB), ALLOCATABLE, DIMENSION(:) :: DX,RDX,RDXN,X_S,DX_WGT,MF_FRAC,PARTICLE_INSERT_CLOCK REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_0 diff --git a/Source/wall.f90 b/Source/wall.f90 index 07f47234348..63bb4d93d81 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -19,7 +19,6 @@ MODULE WALL_ROUTINES !> \brief Main control routine for applying boundary conditions. -!> !> \param T Current time (s) !> \param DT Current time step (s) !> \param NM Mesh number @@ -36,6 +35,7 @@ SUBROUTINE WALL_BC(T,DT,NM) REAL(EB) :: DT_BC,SLIP_COEF INTEGER :: IW,IP,ICF,ITW TYPE(WALL_TYPE), POINTER :: WC +TYPE(THIN_WALL_TYPE), POINTER :: TW TYPE(SURFACE_TYPE), POINTER :: SF TYPE(CFACE_TYPE), POINTER :: CFA TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP @@ -84,9 +84,47 @@ SUBROUTINE WALL_BC(T,DT,NM) ENDIF ENDIF -!$OMP PARALLEL PRIVATE(IW,WC,SF,BC,B1,B2,LP,LPC,CFA,IP,ICF,SLIP_COEF) +! For OpenMP, spin up threads so that they run through all the routines + +!$OMP PARALLEL PRIVATE(IW,WC,TW,SF,BC,B1,B2,LP,LPC,CFA,IP,ICF,SLIP_COEF) + +! Run through thermally-thick WALL and THIN_WALL cells and set gas phase near-surface values and heat transfer coefficient. +! This must be done before the SOLID_HEAT_TRANSFER heat conduction and pyrolysis routine to avoid data races with OpenMP. + +!$OMP DO SCHEDULE(DYNAMIC) +WALL_CELL_LOOP_0: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS + WC=>WALL(IW) + BC => BOUNDARY_COORD(WC%BC_INDEX) + B1 => BOUNDARY_PROP1(WC%B1_INDEX) + IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1) + IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP_0 + SF => SURFACE(WC%SURF_INDEX) + CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,WALL_INDEX=IW) + IF (CALL_HT_1D .AND. SF%THERMAL_BC_INDEX==THERMALLY_THICK) & + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,B1%TMP_G-B1%TMP_F,SF,WALL_INDEX_IN=IW) +ENDDO WALL_CELL_LOOP_0 +!$OMP END DO + +IF (N_THIN_WALL_CELLS>0 .AND. CALL_HT_1D) THEN + !$OMP DO SCHEDULE(DYNAMIC) + THIN_WALL_CELL_LOOP_0: DO ITW=1,N_THIN_WALL_CELLS + TW => THIN_WALL(ITW) + BC => BOUNDARY_COORD(TW%BC_INDEX) + B1 => BOUNDARY_PROP1(TW%B1_INDEX) + SF => SURFACE(TW%SURF_INDEX) + CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,TW=TW,THIN_WALL_INDEX=ITW) + IF (TW%WALL_INDEX_M>0) THEN + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,B1%TMP_G-B1%TMP_F,SF,WALL_INDEX_IN=TW%WALL_INDEX_M) + ELSEIF (TW%WALL_INDEX_P>0) THEN + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,B1%TMP_G-B1%TMP_F,SF,WALL_INDEX_IN=TW%WALL_INDEX_P) + ELSE + B1%HEAT_TRANS_COEF = 0._EB + ENDIF + ENDDO THIN_WALL_CELL_LOOP_0 + !$OMP END DO +ENDIF -! Sweep through all WALL cells and apply boundary conditions +! Sweep through all WALL cells and apply thermal, species and density boundary conditions !$OMP DO SCHEDULE(DYNAMIC) WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS @@ -94,13 +132,10 @@ SUBROUTINE WALL_BC(T,DT,NM) WC=>WALL(IW) BC => BOUNDARY_COORD(WC%BC_INDEX) B1 => BOUNDARY_PROP1(WC%B1_INDEX) - IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1) IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP B2 => BOUNDARY_PROP2(WC%B2_INDEX) SF => SURFACE(WC%SURF_INDEX) - CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,WALL_INDEX=IW) - IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW) ELSEIF (CALL_HT_1D) THEN @@ -148,7 +183,7 @@ SUBROUTINE WALL_BC(T,DT,NM) !$OMP END DO ENDIF -! Loop through all CFACEs and apply heat transfer BCs +! Loop through all CFACEs and apply thermal, species and density BCs !$OMP DO SCHEDULE(DYNAMIC) CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS @@ -162,6 +197,7 @@ SUBROUTINE WALL_BC(T,DT,NM) IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,CFACE_INDEX=ICF) ELSEIF (CALL_HT_1D) THEN + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,B1%TMP_G-B1%TMP_F,SF,CFACE_INDEX_IN=ICF) CALL SOLID_HEAT_TRANSFER(NM,T,DT_BC,CFACE_INDEX=ICF) ENDIF @@ -211,6 +247,7 @@ SUBROUTINE WALL_BC(T,DT,NM) IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,PARTICLE_INDEX=IP) ELSEIF (CALL_HT_1D) THEN + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,B1%TMP_G-B1%TMP_F,SF,PARTICLE_INDEX_IN=IP) CALL SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX=IP) ENDIF @@ -300,15 +337,16 @@ END SUBROUTINE ASSIGN_GHOST_VALUE !> \param WALL_INDEX Index of wall cell !> \param PARTICLE_INDEX Index of particle -SUBROUTINE NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP,WALL_INDEX,PARTICLE_INDEX) +SUBROUTINE NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP,TW,WALL_INDEX,PARTICLE_INDEX,THIN_WALL_INDEX) USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP REAL(EB), INTENT(IN) :: T -INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,PARTICLE_INDEX -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 +INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,PARTICLE_INDEX,THIN_WALL_INDEX +TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1,B1M,B1P TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(SURFACE_TYPE), POINTER :: SF TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER, OPTIONAL :: LP +TYPE(THIN_WALL_TYPE), POINTER, OPTIONAL :: TW REAL(EB) :: TSI,RAMP_FACTOR,UBAR,VBAR,WBAR IF (PRESENT(WALL_INDEX)) THEN @@ -339,13 +377,41 @@ SUBROUTINE NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP,WALL_INDEX,PARTICLE_INDEX) B1%U_TANG = SQRT(UBAR**2+VBAR**2+WBAR**2) ENDIF -IF (SF%TMP_GAS_FRONT > 0._EB) THEN - B1%TMP_G = TMPA + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGF)%INDEX)*(SF%TMP_GAS_FRONT-TMPA) -ELSE - B1%TMP_G = TMP(BC%IIG,BC%JJG,BC%KKG) +! Set near-wall gas temperature, B1%TMP_G, and incoming radiation, B1%Q_RAD_IN + +IF (PRESENT(WALL_INDEX) .OR. PRESENT(PARTICLE_INDEX)) THEN + + IF (SF%TMP_GAS_FRONT > 0._EB) THEN + B1%TMP_G = TMPA + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGF)%INDEX)*(SF%TMP_GAS_FRONT-TMPA) + ELSE + B1%TMP_G = TMP(BC%IIG,BC%JJG,BC%KKG) + ENDIF + B1%RHO_G = RHOP(BC%IIG,BC%JJG,BC%KKG) + B1%ZZ_G(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES) + +ELSEIF (PRESENT(THIN_WALL_INDEX)) THEN + + IF (TW%WALL_INDEX_M>0 .OR. TW%WALL_INDEX_P>0) THEN + IF (TW%WALL_INDEX_M>0) THEN + B1M => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_M)%B1_INDEX) + ELSE + B1M => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_P)%B1_INDEX) + ENDIF + IF (TW%WALL_INDEX_P>0) THEN + B1P => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_P)%B1_INDEX) + ELSE + B1P => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_M)%B1_INDEX) + ENDIF + B1%TMP_G = 0.5_EB*(B1M%TMP_G+B1P%TMP_G) + B1%Q_RAD_IN = 0.5_EB*(B1M%Q_RAD_IN+B1P%Q_RAD_IN) + ELSE + B1%TMP_G = TMP(BC%IIG,BC%JJG,BC%KKG) + B1%Q_RAD_IN = 0._EB + ENDIF + ! Special case where the gas temperature is fixed by the user + IF (SF%TMP_GAS_FRONT > 0._EB) B1%TMP_G = TMPA + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGF)%INDEX)*(SF%TMP_GAS_FRONT-TMPA) + ENDIF -B1%RHO_G = RHOP(BC%IIG,BC%JJG,BC%KKG) -B1%ZZ_G(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES) END SUBROUTINE NEAR_SURFACE_GAS_VARIABLES @@ -375,7 +441,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES),ZZ_GET_OTHER(1:N_TRACKED_SPECIES),& RHO_OTHER,RHO_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),RHO_ZZ_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2, & DDO,PBAR_G,PBAR_OTHER,DENOM,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER, & - MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER,SF_HTC,& + MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER,& BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G,MW_G,MW_G_2,MW_OTHER,MW_OTHER_2 LOGICAL :: SECOND_ORDER_INTERPOLATED_BOUNDARY,ATMOSPHERIC_INTERPOLATION,CC_SOLID_FLAG INTEGER :: IIO,JJO,KKO,N,ADCOUNT,ICG,ICO @@ -412,16 +478,8 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I PY => PROPERTY(LP%PROP_INDEX) IF (PY%HEAT_TRANSFER_COEFFICIENT>0._EB) THEN ! the user has added a PROP line with a specified HTC B1%HEAT_TRANS_COEF = PY%HEAT_TRANSFER_COEFFICIENT - ELSEIF (SF%H_FIXED>=0._EB) THEN ! the user has assigned a SURF_ID to the PART line with HTC>0 - SF_HTC = SF%H_FIXED - IF (SF%RAMP_H_FIXED_INDEX>0) SF_HTC = SF_HTC * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX) - B1%HEAT_TRANS_COEF = SF_HTC - ELSEIF (SF%ID=='MASSLESS TARGET') THEN ! the particle for the gas phase device is not explicitly defined and needs an HTC - B1%HEAT_TRANS_COEF = 10._EB ELSE ! the HTC will be computed using the surrounding gas phase environment - SF_HTC = -1._EB - B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,B1%TMP_G-PY%GAUGE_TEMPERATURE,SF_HTC,SF,& - PARTICLE_INDEX_IN=PARTICLE_INDEX) + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,B1%TMP_G-PY%GAUGE_TEMPERATURE,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX) ENDIF B1%Q_CON_F = B1%HEAT_TRANS_COEF*(B1%TMP_G-PY%GAUGE_TEMPERATURE) RETURN @@ -504,19 +562,12 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I DTMP = B1%TMP_G - B1%TMP_F - IF (SF%H_FIXED >= 0._EB) THEN - SF_HTC = SF%H_FIXED - IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC = SF_HTC * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX) - ELSE - SF_HTC = -1._EB - ENDIF - IF (PRESENT(WALL_INDEX)) THEN - B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,WALL_INDEX_IN=WALL_INDEX) + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,DTMP,SF,WALL_INDEX_IN=WALL_INDEX) ELSEIF (PRESENT(CFACE_INDEX)) THEN - B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,CFACE_INDEX_IN=CFACE_INDEX) + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,DTMP,SF,CFACE_INDEX_IN=CFACE_INDEX) ELSEIF (PRESENT(PARTICLE_INDEX)) THEN - B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX) + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T,DTMP,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX) ENDIF B1%Q_CON_F = B1%HEAT_TRANS_COEF*DTMP IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4 @@ -540,13 +591,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I QNET = -RAMP_FACTOR*SF%CONVECTIVE_HEAT_FLUX*B1%AREA_ADJUST ENDIF - IF (SF%H_FIXED >= 0._EB) THEN - SF_HTC = SF%H_FIXED - IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC = SF_HTC * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX) - ELSE - SF_HTC = -1._EB - ENDIF - IF (ABS(SF_HTC) < TWO_EPSILON_EB) THEN + IF (ABS(SF%H_FIXED) < TWO_EPSILON_EB) THEN IF (RADIATION) B1%TMP_F = ((-QNET + B1%Q_RAD_IN)/(B1%EMISSIVITY * SIGMA))**0.25_EB B1%HEAT_TRANS_COEF = 0._EB B1%Q_CON_F = 0._EB @@ -557,11 +602,11 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I DTMP = B1%TMP_G - B1%TMP_F IF (ABS(QNET) > 0._EB .AND. ABS(DTMP) 0 .OR. TW%WALL_INDEX_P>0) THEN - IF (TW%WALL_INDEX_M>0) THEN - B1M => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_M)%B1_INDEX) - ELSE - B1M => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_P)%B1_INDEX) - ENDIF - IF (TW%WALL_INDEX_P>0) THEN - B1P => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_P)%B1_INDEX) - ELSE - B1P => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_M)%B1_INDEX) - ENDIF - B1%TMP_G = 0.5_EB*(B1M%TMP_G+B1P%TMP_G) - B1%Q_RAD_IN = 0.5_EB*(B1M%Q_RAD_IN+B1P%Q_RAD_IN) - ELSE - B1%TMP_G = TMP(BC%IIG,BC%JJG,BC%KKG) - ISOLATED_THIN_WALL = .TRUE. - ENDIF - ! Special case where the gas temperature is fixed by the user - - IF (SF%TMP_GAS_FRONT > 0._EB) B1%TMP_G = TMPA + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGF)%INDEX)*(SF%TMP_GAS_FRONT-TMPA) + IF (TW%WALL_INDEX_M==0 .AND. TW%WALL_INDEX_P==0) ISOLATED_THIN_WALL = .TRUE. ELSEIF (PRESENT(CFACE_INDEX)) THEN UNPACK_WALL_PARTICLE @@ -1952,22 +1978,6 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, B1%LAYER_REMOVED = .FALSE. -! Evaluate any heat transfer coefficient ramp - -IF (SF%H_FIXED >= 0._EB) THEN - SF_HTC_F = SF%H_FIXED - IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC_F = SF_HTC_F * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX) -ELSE - SF_HTC_F = -1._EB -ENDIF - -IF (SF%H_FIXED_B >= 0._EB) THEN - SF_HTC_B = SF%H_FIXED_B - IF (SF%RAMP_H_FIXED_B_INDEX > 0) SF_HTC_B = SF_HTC_B * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_B_INDEX) -ELSE - SF_HTC_B = -1._EB -ENDIF - ! If the fuel has burned away, return IF (B1%BURNAWAY) THEN @@ -2083,17 +2093,11 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, DTMP = B1%TMP_G - B1%TMP_F IF (PRESENT(WALL_INDEX)) THEN - HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=WALL_INDEX) + HTCF = B1%HEAT_TRANS_COEF ELSEIF (PRESENT(THIN_WALL_INDEX)) THEN - IF (TW%WALL_INDEX_M>0) THEN - HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=TW%WALL_INDEX_M) - ELSEIF (TW%WALL_INDEX_P>0) THEN - HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=TW%WALL_INDEX_P) - ELSE - HTCF = 0._EB - ENDIF + HTCF = B1%HEAT_TRANS_COEF ELSEIF (PRESENT(CFACE_INDEX)) THEN - HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,CFACE_INDEX_IN=CFACE_INDEX) + HTCF = B1%HEAT_TRANS_COEF ELSEIF (PRESENT(PARTICLE_INDEX)) THEN RADIUS = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS)) SELECT CASE(SF%GEOMETRY) @@ -2101,7 +2105,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, CASE (SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(2._EB*DT_BC_SUB) CASE (SURF_SPHERICAL) ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(3._EB*DT_BC_SUB) END SELECT - HTCF = MIN(HTC_LIMIT , HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX)) + HTCF = MIN(HTC_LIMIT , B1%HEAT_TRANS_COEF) ENDIF Q_CON_F = HTCF*DTMP @@ -2125,7 +2129,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, TMP_GAS_BACK = TMP_0(BC%KK) ENDIF DTMP = TMP_GAS_BACK - B1%TMP_B - HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF) + HTCB = HEAT_TRANSFER_COEFFICIENT(NM,T,DTMP,SF,BACK_SIDE=.TRUE.) Q_CON_B = HTCB*DTMP Q_RAD_IN_B = EMISSIVITY_BACK*SIGMA*TMP_GAS_BACK**4 Q_LIQUID_B = 0._EB @@ -2143,34 +2147,13 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, CASE(EXPOSED) ! The backside is exposed to gas in current or adjacent mesh. - Q_LIQUID_B = 0._EB - - IF (BACK_MESH/=NM .AND. BACK_MESH>0) THEN ! Back side is in other mesh. - TMP_GAS_BACK = B1_BACK%TMP_G - DTMP = TMP_GAS_BACK - B1%TMP_B - IF (.NOT.ISOLATED_THIN_WALL_BACK) THEN - HTCB = HEAT_TRANSFER_COEFFICIENT(BACK_MESH,DTMP,SF_HTC_B,SF_BACK,WALL_INDEX_IN=BACK_WALL_INDEX) - ELSE - HTCB = 0._EB - ENDIF - Q_RAD_IN_B = B1_BACK%Q_RAD_IN - ELSE ! Back side is in current mesh. - TMP_GAS_BACK = TMP(BC_BACK%IIG,BC_BACK%JJG,BC_BACK%KKG) - DTMP = TMP_GAS_BACK - B1%TMP_B - IF (PRESENT(WALL_INDEX)) & - HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF_BACK,WALL_INDEX_IN=ONE_D%BACK_INDEX) - IF (PRESENT(THIN_WALL_INDEX)) THEN - IF (BACK_WALL_INDEX>0) THEN - HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF_BACK,WALL_INDEX_IN=BACK_WALL_INDEX) - ELSE - HTCB = 0._EB - ENDIF - ENDIF - IF (PRESENT(CFACE_INDEX)) & - HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF_BACK,CFACE_INDEX_IN=ONE_D%BACK_INDEX) - B1_BACK%HEAT_TRANS_COEF = HTCB - Q_RAD_IN_B = B1_BACK%Q_RAD_IN - IF (N_LP_ARRAY_INDICES>0 .AND. PRESENT(WALL_INDEX)) Q_LIQUID_B = -SUM(B2_BACK%LP_CPUA(:)) + B1%Q_CONDENSE + TMP_GAS_BACK = B1_BACK%TMP_G + HTCB = B1_BACK%HEAT_TRANS_COEF + Q_RAD_IN_B = B1_BACK%Q_RAD_IN + IF (N_LP_ARRAY_INDICES>0 .AND. PRESENT(WALL_INDEX)) THEN + Q_LIQUID_B = -SUM(B2_BACK%LP_CPUA(:)) + B1%Q_CONDENSE + ELSE + Q_LIQUID_B = 0._EB ENDIF Q_CON_B = HTCB*DTMP @@ -2231,9 +2214,9 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, IF (SF%SURFACE_VOLUME_RATIO(LAYER_INDEX(I))<=0._EB) CYCLE DTMP = B1%TMP_G - ONE_D%TMP(I) IF (PRESENT(WALL_INDEX)) THEN - HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=WALL_INDEX) + HTCF = HEAT_TRANSFER_COEFFICIENT(NM,T,DTMP,SF,WALL_INDEX_IN=WALL_INDEX) ELSEIF(PRESENT(CFACE_INDEX)) THEN - HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,CFACE_INDEX_IN=CFACE_INDEX) + HTCF = HEAT_TRANSFER_COEFFICIENT(NM,T,DTMP,SF,CFACE_INDEX_IN=CFACE_INDEX) ENDIF DEL_DOT_Q_SC = HTCF*DTMP Q_ADD(I) = Q_ADD(I) + SF%SURFACE_VOLUME_RATIO(LAYER_INDEX(I))*SF%PACKING_RATIO(LAYER_INDEX(I))*DEL_DOT_Q_SC @@ -2860,7 +2843,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, B1%Q_CON_F = B1%Q_CON_F / DT_BC IF (RADIATION .AND. .NOT. ONE_D%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT / DT_BC * SIGMA * B1%EMISSIVITY -IF (.NOT. SF%BOUNDARY_FUEL_MODEL) B1%HEAT_TRANS_COEF = HTCF +! IF (.NOT. SF%BOUNDARY_FUEL_MODEL) B1%HEAT_TRANS_COEF = HTCF ! If any gas massflux or particle mass flux is non-zero or the surface temperature exceeds the ignition temperature, ! set the ignition time @@ -3112,7 +3095,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F REAL(EB) :: REACTION_RATE,Y_O2,X_O2,MW(N_MATS),Y_GAS(N_MATS),Y_TMP(N_MATS),Y_SV(N_MATS),X_SV(N_MATS),X_L(N_MATS),& D_FILM,H_MASS,RE_L,SHERWOOD,MFLUX,MU_FILM,SC_FILM,TMP_FILM,TMP_G,U2,V2,W2,VEL,& DR,R_S_0,R_S_1,H_R,H_R_B,H_S_B,H_S,LENGTH_SCALE,SUM_Y_GAS,SUM_Y_SV,NU_O2_CHAR,Y_O2_S,& - SUM_Y_SV_SMIX(N_TRACKED_SPECIES),X_L_SUM,RHO_DOT_EXTRA,MFLUX_MAX,RHO_FILM,CP_FILM,PR_FILM,K_FILM,EVAP_FILM_FAC,& + SUM_Y_SV_SMIX(N_TRACKED_SPECIES),X_L_SUM,RHO_DOT_EXTRA,MFLUX_MAX,RHO_FILM,CP_FILM,PR_FILM,K_FILM,& RHO_DOT,RHO_DOT_REAC(MAX_REACTIONS),RHO_DOT_REAC_SUM,H_MASS_DNS LOGICAL :: LIQUID(N_MATS),SPEC_ID_ALREADY_USED(N_MATS),DO_EVAPORATION @@ -3155,6 +3138,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F IF (X_L_SUM < TWO_EPSILON_EB) RETURN + Y_GAS = 0._EB SUM_Y_GAS = 0._EB SPEC_ID_ALREADY_USED = .FALSE. SMIX_INDEX = 0 @@ -3226,16 +3210,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F ENDDO MATERIAL_LOOP_2 Y_GAS = Y_TMP - ! Get film properties - - SELECT CASE(SF%GEOMETRY) - CASE DEFAULT ; EVAP_FILM_FAC = PLATE_FILM_FAC - CASE(SURF_SPHERICAL) ; EVAP_FILM_FAC = SPHERE_FILM_FAC - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; EVAP_FILM_FAC = PLATE_FILM_FAC - CASE(SURF_CARTESIAN) ; EVAP_FILM_FAC = PLATE_FILM_FAC - END SELECT - - CALL GET_FILM_PROPERTIES(N_MATS,EVAP_FILM_FAC,Y_SV,Y_GAS,SMIX_INDEX,TMP_F,TMP(IIG,JJG,KKG),ZZ_GET,& + CALL GET_FILM_PROPERTIES(N_MATS,SF%FILM_FACTOR,Y_SV,Y_GAS,SMIX_INDEX,TMP_F,TMP(IIG,JJG,KKG),ZZ_GET,& PBAR(KKG,PRESSURE_ZONE(IIG,JJG,KKG)),TMP_FILM,MU_FILM,K_FILM,CP_FILM,D_FILM,& RHO_FILM,PR_FILM,SC_FILM) @@ -3500,67 +3475,70 @@ END SUBROUTINE PYROLYSIS !> \brief Compute the convective heat transfer coefficient +!> \param NM Mesh number +!> \param T Current time (s) +!> \param DELTA_N_TMP Difference between gas and surface temperature (K) +!> \param SF SURFACE derived type variable +!> \param WALL_INDEX_IN Optional wall cell index +!> \param CFACE_INDEX_IN Optional cface index +!> \param PARTICLE_INDEX_IN Optional particle index +!> \param BACK_SIZE Optional flag indicating if the surface is on the back side of the obstruction -REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_INDEX_IN,CFACE_INDEX_IN,PARTICLE_INDEX_IN) +REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NM,T,DELTA_N_TMP,SF,WALL_INDEX_IN,CFACE_INDEX_IN,PARTICLE_INDEX_IN,BACK_SIDE) USE TURBULENCE, ONLY: LOGLAW_HEAT_FLUX_MODEL,NATURAL_CONVECTION_MODEL,FORCED_CONVECTION_MODEL,RAYLEIGH_HEAT_FLUX_MODEL,& FM_HEAT_FLUX_MODEL USE PHYSICAL_FUNCTIONS, ONLY: GET_CONDUCTIVITY,GET_VISCOSITY,GET_SPECIFIC_HEAT -REAL(EB), INTENT(IN) :: DELTA_N_TMP,H_FIXED -INTEGER, INTENT(IN) :: NMX +REAL(EB), INTENT(IN) :: DELTA_N_TMP,T +INTEGER, INTENT(IN) :: NM INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX_IN,PARTICLE_INDEX_IN,CFACE_INDEX_IN +LOGICAL, INTENT(IN), OPTIONAL :: BACK_SIDE 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,XX + R_DROP,CONV_LENGTH,GR,RA,NUSSELT_FORCED,NUSSELT_FREE,NUSSELT_IMPINGE,PHI,XX,H_FIXED INTEGER, PARAMETER :: NATURAL=1,FORCED=2,IMPACT=3,RESOLVED=4 -TYPE(MESH_TYPE), POINTER :: MX -TYPE(SURFACE_TYPE), INTENT(IN), POINTER :: SFX -TYPE(WALL_TYPE), POINTER :: WCX -TYPE(CFACE_TYPE), POINTER :: CFAX -TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LPX -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: P1X -TYPE(BOUNDARY_PROP2_TYPE), POINTER :: P2X -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BCX - -MX => MESHES(NMX) -CONV_LENGTH = SFX%CONV_LENGTH +TYPE(MESH_TYPE), POINTER :: M +TYPE(SURFACE_TYPE), INTENT(IN), POINTER :: SF +TYPE(WALL_TYPE), POINTER :: WC +TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP +TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 +TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2 +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC + +M => MESHES(NM) +CONV_LENGTH = SF%CONV_LENGTH +IF (PRESENT(BACK_SIDE)) THEN + H_FIXED = SF%H_FIXED_B +ELSE + H_FIXED = SF%H_FIXED +ENDIF ! Determine if this is a particle or wall cell IF (PRESENT(PARTICLE_INDEX_IN)) THEN - LPX => LAGRANGIAN_PARTICLE(PARTICLE_INDEX_IN) - P1X => BOUNDARY_PROP1(LPX%B1_INDEX) - IF (LAGRANGIAN_PARTICLE_CLASS(LPX%CLASS_INDEX)%INCLUDE_BOUNDARY_PROP2_TYPE) & - P2X => BOUNDARY_PROP2(LPX%B2_INDEX) - BCX => BOUNDARY_COORD(LPX%BC_INDEX) - DN = SFX%CONV_LENGTH - R_DROP = LPX%RADIUS - IF (R_DROP>TWO_EPSILON_EB .AND. SFX%GEOMETRY/=SURF_CARTESIAN) CONV_LENGTH = 2._EB*R_DROP - TMP_G = P1X%TMP_G - RHO_G = P1X%RHO_G - ZZ_G(1:N_TRACKED_SPECIES) = P1X%ZZ_G(1:N_TRACKED_SPECIES) + LP => M%LAGRANGIAN_PARTICLE(PARTICLE_INDEX_IN) + B1 => M%BOUNDARY_PROP1(LP%B1_INDEX) + IF (LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)%INCLUDE_BOUNDARY_PROP2_TYPE) B2 => M%BOUNDARY_PROP2(LP%B2_INDEX) + BC => M%BOUNDARY_COORD(LP%BC_INDEX) + DN = SF%CONV_LENGTH + R_DROP = LP%RADIUS + IF (R_DROP>TWO_EPSILON_EB .AND. SF%GEOMETRY/=SURF_CARTESIAN) CONV_LENGTH = 2._EB*R_DROP ELSEIF (PRESENT(WALL_INDEX_IN)) THEN - WCX => MX%WALL(WALL_INDEX_IN) - P1X => MX%BOUNDARY_PROP1(WCX%B1_INDEX) - P2X => MX%BOUNDARY_PROP2(WCX%B2_INDEX) - BCX => MX%BOUNDARY_COORD(WCX%BC_INDEX) - DN = 1._EB/P1X%RDN - TMP_G = P1X%TMP_G - RHO_G = P1X%RHO_G - ZZ_G(1:N_TRACKED_SPECIES) = P1X%ZZ_G(1:N_TRACKED_SPECIES) + WC => M%WALL(WALL_INDEX_IN) + B1 => M%BOUNDARY_PROP1(WC%B1_INDEX) + B2 => M%BOUNDARY_PROP2(WC%B2_INDEX) + BC => M%BOUNDARY_COORD(WC%BC_INDEX) + DN = 1._EB/B1%RDN ELSEIF (PRESENT(CFACE_INDEX_IN)) THEN - CFAX => CFACE(CFACE_INDEX_IN) - P1X => BOUNDARY_PROP1(CFAX%B1_INDEX) - P2X => BOUNDARY_PROP2(CFAX%B2_INDEX) - BCX => BOUNDARY_COORD(CFAX%BC_INDEX) - DN = 1._EB/P1X%RDN - TMP_G = P1X%TMP_G - RHO_G = P1X%RHO_G - ZZ_G(1:N_TRACKED_SPECIES) = P1X%ZZ_G(1:N_TRACKED_SPECIES) + CFA => M%CFACE(CFACE_INDEX_IN) + B1 => M%BOUNDARY_PROP1(CFA%B1_INDEX) + B2 => M%BOUNDARY_PROP2(CFA%B2_INDEX) + BC => M%BOUNDARY_COORD(CFA%BC_INDEX) + DN = 1._EB/B1%RDN ELSE IF (H_FIXED >= 0._EB) THEN - HEAT_TRANSFER_COEFFICIENT = H_FIXED + CALL CONSTANT_HTC ELSE HEAT_TRANSFER_COEFFICIENT = 1.31_EB*ABS(DELTA_N_TMP)**ONTH ! Natural convection for vertical plane, Holman, 10th, Tab. 7.2 ENDIF @@ -3569,20 +3547,13 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND ! Calculate HEAT_TRANSFER_COEFFICIENT -SELECT CASE(SFX%GEOMETRY) - CASE DEFAULT ; FILM_FAC = PLATE_FILM_FAC - CASE(SURF_SPHERICAL) ; FILM_FAC = SPHERE_FILM_FAC - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; FILM_FAC = PLATE_FILM_FAC - CASE(SURF_CARTESIAN) ; FILM_FAC = PLATE_FILM_FAC -END SELECT - -TMP_FILM = P1X%TMP_F+FILM_FAC*(TMP_G-P1X%TMP_F) +TMP_FILM = B1%TMP_F+SF%FILM_FACTOR*(B1%TMP_G-B1%TMP_F) ! If the user wants a specified HTC, set it and return -H_FIXED_IF: IF (H_FIXED >= 0._EB .AND. SFX%HEAT_TRANSFER_MODEL/=IMPINGING_JET_HTC_MODEL) THEN +H_FIXED_IF: IF (H_FIXED >= 0._EB .AND. SF%HEAT_TRANSFER_MODEL/=IMPINGING_JET_HTC_MODEL) THEN - HEAT_TRANSFER_COEFFICIENT = H_FIXED + CALL CONSTANT_HTC ELSE H_FIXED_IF @@ -3591,82 +3562,82 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND HTC_TYPE_IF: IF ( (SIM_MODE==DNS_MODE .OR. SOLID_PHASE_ONLY) .AND. (PRESENT(WALL_INDEX_IN) .OR. PRESENT(CFACE_INDEX_IN)) ) THEN - HEAT_TRANSFER_COEFFICIENT = 2._EB * P1X%K_G * P1X%RDN + HEAT_TRANSFER_COEFFICIENT = 2._EB * B1%K_G * B1%RDN IF (SIM_MODE==DNS_MODE) RETURN ELSE HTC_TYPE_IF - CALL GET_VISCOSITY(ZZ_G,MU_G,TMP_FILM) - CALL GET_CONDUCTIVITY(ZZ_G,K_G,TMP_FILM) + CALL GET_VISCOSITY(B1%ZZ_G,MU_G,TMP_FILM) + CALL GET_CONDUCTIVITY(B1%ZZ_G,K_G,TMP_FILM) - HTC_MODEL_SELECT: SELECT CASE(SFX%HEAT_TRANSFER_MODEL) + HTC_MODEL_SELECT: SELECT CASE(SF%HEAT_TRANSFER_MODEL) CASE(DEFAULT_HTC_MODEL,IMPINGING_JET_HTC_MODEL) - RE = RHO_G*P1X%U_TANG*CONV_LENGTH/MU_G - GR = GRAV*ABS(DELTA_N_TMP)*CONV_LENGTH**3*(RHO_G/MU_G)**2/TMP_FILM - IF (SFX%BOUNDARY_FUEL_MODEL) THEN + RE = B1%RHO_G*B1%U_TANG*CONV_LENGTH/MU_G + GR = GRAV*ABS(DELTA_N_TMP)*CONV_LENGTH**3*(B1%RHO_G/MU_G)**2/TMP_FILM + IF (SF%BOUNDARY_FUEL_MODEL) THEN SURF_GEOMETRY = SURF_CYLINDRICAL ELSE - SURF_GEOMETRY = SFX%GEOMETRY + SURF_GEOMETRY = SF%GEOMETRY ENDIF ! Check if custom Nusselt correlation is defined - IF (ANY((/SFX%NUSSELT_C0,SFX%NUSSELT_C1,SFX%NUSSELT_C2,SFX%NUSSELT_M/)>0._EB)) THEN + IF (ANY((/SF%NUSSELT_C0,SF%NUSSELT_C1,SF%NUSSELT_C2,SF%NUSSELT_M/)>0._EB)) THEN CALL FORCED_CONVECTION_MODEL(NUSSELT_FORCED,RE,PR_ONTH,SURF_GEOMETRY,& - SFX%NUSSELT_C0,SFX%NUSSELT_C1,SFX%NUSSELT_C2,SFX%NUSSELT_M) + SF%NUSSELT_C0,SF%NUSSELT_C1,SF%NUSSELT_C2,SF%NUSSELT_M) ELSE CALL FORCED_CONVECTION_MODEL(NUSSELT_FORCED,RE,PR_ONTH,SURF_GEOMETRY) ENDIF RA = GR*PR_AIR - CALL NATURAL_CONVECTION_MODEL(NUSSELT_FREE,RA,SFX,BCX%IOR,DELTA_N_TMP) + CALL NATURAL_CONVECTION_MODEL(NUSSELT_FREE,RA,SF,BC%IOR,DELTA_N_TMP) NUSSELT_IMPINGE = 0._EB - IF (SFX%HEAT_TRANSFER_MODEL==IMPINGING_JET_HTC_MODEL) THEN - 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 + IF (SF%HEAT_TRANSFER_MODEL==IMPINGING_JET_HTC_MODEL) THEN + XX = SQRT( (SF%XYZ(1)-BC%X)**2 + (SF%XYZ(2)-BC%Y)**2 + (SF%XYZ(3)-BC%Z)**2 ) + NUSSELT_IMPINGE = HEAT_TRANSFER_COEFFICIENT*EXP(-0.5_EB * (XX / SF%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 ! BOUNDARY_PROP2 not allocated for particles ELSE - P2X%HEAT_TRANSFER_REGIME=MAXLOC((/NUSSELT_FREE,NUSSELT_FORCED,NUSSELT_IMPINGE,2._EB*CONV_LENGTH/DN/),DIM=1) - SELECT CASE(P2X%HEAT_TRANSFER_REGIME) + B2%HEAT_TRANSFER_REGIME=MAXLOC((/NUSSELT_FREE,NUSSELT_FORCED,NUSSELT_IMPINGE,2._EB*CONV_LENGTH/DN/),DIM=1) + SELECT CASE(B2%HEAT_TRANSFER_REGIME) CASE(NATURAL) - P2X%Z_STAR = NUSSELT_FREE**0.25_EB + B2%Z_STAR = NUSSELT_FREE**0.25_EB CASE(FORCED) - P2X%Z_STAR = SQRT(NUSSELT_FORCED) + B2%Z_STAR = SQRT(NUSSELT_FORCED) CASE(IMPACT) - P2X%Z_STAR = SQRT(NUSSELT_IMPINGE) + B2%Z_STAR = SQRT(NUSSELT_IMPINGE) CASE(RESOLVED) - P2X%Z_STAR = 1._EB + B2%Z_STAR = 1._EB END SELECT HEAT_TRANSFER_COEFFICIENT = MAX(NUSSELT_FREE,NUSSELT_FORCED,NUSSELT_IMPINGE,2._EB*CONV_LENGTH/DN)*K_G/CONV_LENGTH ENDIF CASE(LOGLAW_HTC_MODEL) - CALL GET_SPECIFIC_HEAT(ZZ_G,CP_G,TMP_FILM) - FRICTION_VELOCITY = P2X%U_TAU - YPLUS = P2X%Y_PLUS - CALL LOGLAW_HEAT_FLUX_MODEL(H_FORCED,YPLUS,FRICTION_VELOCITY,K_G,RHO_G,CP_G,MU_G) + CALL GET_SPECIFIC_HEAT(B1%ZZ_G,CP_G,TMP_FILM) + FRICTION_VELOCITY = B2%U_TAU + YPLUS = B2%Y_PLUS + CALL LOGLAW_HEAT_FLUX_MODEL(H_FORCED,YPLUS,FRICTION_VELOCITY,K_G,B1%RHO_G,CP_G,MU_G) HEAT_TRANSFER_COEFFICIENT = H_FORCED 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) - P2X%Z_STAR = ZSTAR - P2X%HEAT_TRANSFER_REGIME = HTR + CALL GET_SPECIFIC_HEAT(B1%ZZ_G,CP_G,TMP_FILM) + CALL RAYLEIGH_HEAT_FLUX_MODEL(H_NATURAL,ZSTAR,HTR,DN,B1%TMP_F,B1%TMP_G,K_G,B1%RHO_G,CP_G,MU_G,B1%U_TANG) + B2%Z_STAR = ZSTAR + B2%HEAT_TRANSFER_REGIME = HTR HEAT_TRANSFER_COEFFICIENT = H_NATURAL CASE(FM_HTC_MODEL) - CALL GET_SPECIFIC_HEAT(ZZ_G,CP_G,TMP_FILM) - CALL FM_HEAT_FLUX_MODEL(H_NATURAL,DN,P1X%TMP_F,TMP_G,K_G,RHO_G,CP_G) + CALL GET_SPECIFIC_HEAT(B1%ZZ_G,CP_G,TMP_FILM) + CALL FM_HEAT_FLUX_MODEL(H_NATURAL,DN,B1%TMP_F,B1%TMP_G,K_G,B1%RHO_G,CP_G) HEAT_TRANSFER_COEFFICIENT = H_NATURAL END SELECT HTC_MODEL_SELECT ENDIF HTC_TYPE_IF ENDIF H_FIXED_IF -IF (SFX%BLOWING .AND. .NOT. SFX%BOUNDARY_FUEL_MODEL .AND. SIM_MODE /= DNS_MODE .AND. ALLOCATED(P1X%M_DOT_G_PP_ACTUAL)) THEN +IF (SF%BLOWING .AND. .NOT. SF%BOUNDARY_FUEL_MODEL .AND. SIM_MODE /= DNS_MODE .AND. ALLOCATED(B1%M_DOT_G_PP_ACTUAL)) THEN PHI = 0._EB - IF (SFX%INCLUDE_BOUNDARY_PROP2_TYPE) P2X%BLOWING_CORRECTION=0._EB + IF (SF%INCLUDE_BOUNDARY_PROP2_TYPE) B2%BLOWING_CORRECTION=0._EB ITMP = INT(TMP_FILM) DO I=1,N_TRACKED_SPECIES - IF (ABS(P1X%M_DOT_G_PP_ACTUAL(I)) <= TWO_EPSILON_EB) CYCLE - PHI = PHI + P1X%M_DOT_G_PP_ACTUAL(I)*CP_Z(ITMP,I) + IF (ABS(B1%M_DOT_G_PP_ACTUAL(I)) <= TWO_EPSILON_EB) CYCLE + PHI = PHI + B1%M_DOT_G_PP_ACTUAL(I)*CP_Z(ITMP,I) ENDDO IF (ABS(PHI)>TWO_EPSILON_EB .AND. ABS(HEAT_TRANSFER_COEFFICIENT)>=TWO_EPSILON_EB) THEN PHI = PHI / HEAT_TRANSFER_COEFFICIENT @@ -3674,11 +3645,26 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND HEAT_TRANSFER_COEFFICIENT = 0._EB ELSE HEAT_TRANSFER_COEFFICIENT = HEAT_TRANSFER_COEFFICIENT * PHI/(EXP(PHI)-1._EB) - IF (SFX%INCLUDE_BOUNDARY_PROP2_TYPE) P2X%BLOWING_CORRECTION = PHI/(EXP(PHI)-1._EB) + IF (SF%INCLUDE_BOUNDARY_PROP2_TYPE) B2%BLOWING_CORRECTION = PHI/(EXP(PHI)-1._EB) ENDIF ENDIF ENDIF +CONTAINS + +SUBROUTINE CONSTANT_HTC + +USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP + +HEAT_TRANSFER_COEFFICIENT = SF%H_FIXED +IF (SF%RAMP_H_FIXED_INDEX>0) THEN + HEAT_TRANSFER_COEFFICIENT = HEAT_TRANSFER_COEFFICIENT*EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX) +ELSE + HEAT_TRANSFER_COEFFICIENT = SF%H_FIXED +ENDIF + +END SUBROUTINE CONSTANT_HTC + END FUNCTION HEAT_TRANSFER_COEFFICIENT @@ -3744,7 +3730,9 @@ SUBROUTINE HT3D_TEMPERATURE_EXCHANGE(NM) ONE_D%TMP(I) = (ONE_D%RHO_C_S(I)*ONE_D%TMP(I) + & THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%RHO_C_S(I_NODE)*ONE_D2%DELTA_TMP(I_NODE))/ONE_D%RHO_C_S(I) ENDDO WEIGHT_LOOP + !$OMP CRITICAL IF (THR_D%NODE(I)%MESH_NUMBER==NM) M%TMP(THR_D%NODE(I)%I,THR_D%NODE(I)%J,THR_D%NODE(I)%K) = ONE_D%TMP(I) + !$OMP END CRITICAL ENDDO NODE_LOOP ONE_D%TMP(0) = ONE_D%TMP(0) + ONE_D%TMP(1) - TMP_1 ONE_D%TMP(NWP+1) = ONE_D%TMP(NWP+1) + ONE_D%TMP(NWP) - TMP_NWP @@ -3855,8 +3843,10 @@ SUBROUTINE TGA_ANALYSIS(NM) T_TGA = REAL(I,EB)*TGA_DT B1%TMP_G = TMPA + TGA_HEATING_RATE*T_TGA IF (TGA_WALL_INDEX>0) THEN + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T_TGA,B1%TMP_G-B1%TMP_F,SF,WALL_INDEX_IN=IW) CALL SOLID_HEAT_TRANSFER(NM,T_TGA,TGA_DT,WALL_INDEX=IW) ELSE + B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,T_TGA,B1%TMP_G-B1%TMP_F,SF,PARTICLE_INDEX_IN=IP) CALL SOLID_HEAT_TRANSFER(NM,T_TGA,TGA_DT,PARTICLE_INDEX=IP) ENDIF IF (I==1 .OR. B1%TMP_G >= TMP_DUMP) THEN diff --git a/Verification/Heat_Transfer/insulated_steel_pipe_2d.fds b/Verification/Heat_Transfer/insulated_steel_pipe_2d.fds index 98941a4e28e..5d908e9d7a1 100644 --- a/Verification/Heat_Transfer/insulated_steel_pipe_2d.fds +++ b/Verification/Heat_Transfer/insulated_steel_pipe_2d.fds @@ -26,7 +26,6 @@ TMP_GAS_FRONT = 480. TMP_GAS_BACK = 20. HEAT_TRANSFER_COEFFICIENT = 10. - HEAT_TRANSFER_COEFFICIENT_BACK = 10. GEOMETRY = 'CYLINDRICAL' INNER_RADIUS = 0.01 THICKNESS = 0.02,0.01,0.02 / @@ -36,7 +35,6 @@ TMP_GAS_FRONT = 20. TMP_GAS_BACK = 480. HEAT_TRANSFER_COEFFICIENT = 10. - HEAT_TRANSFER_COEFFICIENT_BACK = 10. GEOMETRY = 'INNER CYLINDRICAL' INNER_RADIUS = 0.01 THICKNESS = 0.02,0.01,0.02 /