Merge pull request #14271 from mcgratta/master

FDS Source: Issue #14269. Radiation emission output

Author: mcgratta

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -11405,7 +11405,9 @@ \section{Gas Phase Output Quantities}
 \ct{PRESSURE ZONE}                             & Section~\ref{info:ZONE}                           &                   & D,S          \\ \hline
 \ct{Q CRITERION}                               & $\frac{1}{2} [\mathrm{trace}(\nabla \mathbf{u})^2 - \mathrm{trace}( (\nabla \mathbf{u})^2)]$ & 1/s$^2$ & D,I,P,S \\ \hline
 \ct{RADIAL VELOCITY}                           & $(u,v)\cdot(x,y)/\sqrt{x^2+y^2}$                  & m/s               & D,I,P,S      \\ \hline
-\ct{RADIATION LOSS}                            & $\nabla \cdot \bq_r''$                            & kW/m$^3$          & D,I,P,S      \\ \hline
+\ct{RADIATION ABSORPTION}                      & $\kappa U$                                        & kW/m$^3$          & D,I,P,S      \\ \hline
+\ct{RADIATION EMISSION}                        & $4 \kappa \sigma T^4$                             & kW/m$^3$          & D,I,P,S      \\ \hline
+\ct{RADIATION LOSS}                            & $\nabla \cdot \bq_r''=\kappa (U-4 \sigma T^4)$    & kW/m$^3$          & D,I,P,S      \\ \hline
 \ct{RADIATIVE HEAT FLUX GAS}                   & Section~\ref{info:heat_flux}                      & \unit{kW/m^2}          & D            \\ \hline
 \ct{REAC SOURCE TERM}$^1$                      & $\dot{m}_\alpha^{\prime\prime\prime}$             & kg/m$^3$          & D,I,P,S      \\ \hline
 \ct{RELATIVE HUMIDITY}                         & Section~\ref{info:simple_chemistry}               & \%                & D,I,P,S      \\ \hline

Source/cons.f90

@@ -256,6 +256,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: POSITIVE_ERROR_TEST=.FALSE.
 LOGICAL :: OBST_SHAPE_AREA_ADJUST=.FALSE.
 LOGICAL :: STORE_SPECIES_FLUX=.FALSE.
+LOGICAL :: STORE_RADIATION_TERMS=.FALSE.
 LOGICAL :: STORE_PRESSURE_POISSON_RESIDUAL=.FALSE.
 LOGICAL :: OXIDATION_REACTION=.FALSE.
 LOGICAL :: PERIODIC_DOMAIN_X=.FALSE.                !< The domain is periodic \f$ x \f$

Source/data.f90

@@ -439,6 +439,16 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(77)%UNITS = ' '
 OUTPUT_QUANTITY(77)%SHORT_NAME = 'phi'
 
+! Special radiation terms
+
+OUTPUT_QUANTITY(78)%NAME = 'RADIATION EMISSION'
+OUTPUT_QUANTITY(78)%UNITS = 'kW/m3'
+OUTPUT_QUANTITY(78)%SHORT_NAME = 'emit'
+
+OUTPUT_QUANTITY(79)%NAME = 'RADIATION ABSORPTION'
+OUTPUT_QUANTITY(79)%UNITS = 'kW/m3'
+OUTPUT_QUANTITY(79)%SHORT_NAME = 'abs'
+
 ! Cell indices (for Smokeview debugging)
 
 OUTPUT_QUANTITY(80)%NAME = 'CELL INDEX I'

Source/dump.f90

@@ -7618,6 +7618,10 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
 
    CASE(77)  ! LEVEL SET VALUE
       GAS_PHASE_OUTPUT_RES = PHI_LS(II,JJ)
+   CASE(78)  ! RADIATION EMISSION
+      GAS_PHASE_OUTPUT_RES = RADIATION_EMISSION(II,JJ,KK)*0.001_EB
+   CASE(79)  ! RADIATION ABSORPTION
+      GAS_PHASE_OUTPUT_RES = RADIATION_ABSORPTION(II,JJ,KK)*0.001_EB
    CASE(80)  ! CELL INDEX I
       GAS_PHASE_OUTPUT_RES = REAL(II,EB)
    CASE(81)  ! CELL INDEX J

Source/init.f90

@@ -291,6 +291,13 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 ALLOCATE(M%KAPPA_GAS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','KAPPA_GAS',IZERO) ; M%KAPPA_GAS = 0._EB
 ALLOCATE(M%UII(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)       ; CALL ChkMemErr('INIT','UII',IZERO)       ; M%UII       = 0._EB
 
+IF (STORE_RADIATION_TERMS) THEN
+   ALLOCATE(M%RADIATION_EMISSION(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','RADIATION_EMISSION',IZERO)
+   M%RADIATION_EMISSION = 0._EB
+   ALLOCATE(M%RADIATION_ABSORPTION(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','RADIATION_ABSORPTION',IZERO)
+   M%RADIATION_ABSORPTION = 0._EB
+ENDIF
+
 ! Work arrays
 
 ALLOCATE(M%WORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK1',IZERO)

Source/mesh.f90

@@ -46,6 +46,8 @@ MODULE MESH_VARIABLES
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: CHI_R   !< Radiative fraction, \f$ \chi_{{\rm r},ijk} \f$
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: QR      !< Radiation source term, \f$ -\nabla \cdot \dot{\mathbf{q}}_{\rm r}'' \f$
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: QR_W    !< Radiation source term, particles and droplets
+   REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: RADIATION_EMISSION   !< Radiation source term, emission
+   REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: RADIATION_ABSORPTION !< Radiation source term, absorption
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: UII     !< Integrated intensity, \f$ U_{ijk}=\sum_{l=1}^N I_{ijk}^l\delta\Omega^l\f$
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: RSUM    !< \f$ R_0 \sum_\alpha Z_{\alpha,ijk}/W_\alpha \f$
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: D_SOURCE!< Source terms in the expression for the divergence
@@ -348,7 +350,7 @@ MODULE MESH_POINTERS
 
 REAL(EB), POINTER, DIMENSION(:,:,:) :: &
    U,V,W,US,VS,WS,DDDT,D,DS,H,HS,KRES,FVX,FVY,FVZ,FVX_B,FVY_B,FVZ_B,FVX_D,FVY_D,FVZ_D,RHO,RHOS, &
-   MU,MU_DNS,TMP,Q,KAPPA_GAS,CHI_R,QR,QR_W,UII,RSUM,D_SOURCE, &
+   MU,MU_DNS,TMP,Q,KAPPA_GAS,CHI_R,QR,QR_W,RADIATION_EMISSION,RADIATION_ABSORPTION,UII,RSUM,D_SOURCE, &
    CSD2,MTR,MSR,WEM,MIX_TIME,CHEM_SUBIT,STRAIN_RATE,D_Z_MAX,PP_RESIDUAL,LES_FILTER_WIDTH
 REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZ,ZZS,REAC_SOURCE_TERM,DEL_RHO_D_DEL_Z,FX,FY,FZ, &
                                          SWORK1,SWORK2,SWORK3,SWORK4, &
@@ -525,6 +527,8 @@ SUBROUTINE POINT_TO_MESH(NM)
 CHI_R => M%CHI_R
 QR=>M%QR
 QR_W=>M%QR_W
+RADIATION_EMISSION=>M%RADIATION_EMISSION
+RADIATION_ABSORPTION=>M%RADIATION_ABSORPTION
 KAPPA_GAS=>M%KAPPA_GAS
 UII=>M%UII
 PP_RESIDUAL=>M%PP_RESIDUAL

Source/radi.f90

@@ -3546,6 +3546,10 @@ SUBROUTINE RADIATION_FVM
 
 IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
    QR = 0._EB
+   IF (STORE_RADIATION_TERMS) THEN
+      RADIATION_EMISSION   = 0._EB
+      RADIATION_ABSORPTION = 0._EB
+   ENDIF
    IF (N_LP_ARRAY_INDICES>0 .AND. PARTICLES_EXISTED) QR_W = 0._EB
 ENDIF
 
@@ -4476,6 +4480,10 @@ SUBROUTINE RADIATION_FVM
 
    IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
       QR = QR + KAPPA_GAS*UIID(:,:,:,IBND)-KFST4_GAS
+      IF (STORE_RADIATION_TERMS) THEN
+         RADIATION_EMISSION   = RADIATION_EMISSION + KFST4_GAS
+         RADIATION_ABSORPTION = RADIATION_ABSORPTION + KAPPA_GAS*UIID(:,:,:,IBND)
+      ENDIF
       IF (N_LP_ARRAY_INDICES>0 .AND. PARTICLES_EXISTED) QR_W = QR_W + KAPPA_PART*UIID(:,:,:,IBND) - KFST4_PART
    ENDIF
 
@@ -4506,6 +4514,10 @@ SUBROUTINE RADIATION_FVM
 IF (.NOT. (WIDE_BAND_MODEL .OR. WSGG_MODEL)) THEN
    QR = KAPPA_GAS*UII - KFST4_GAS
    IF (N_LP_ARRAY_INDICES>0) QR_W = QR_W + KAPPA_PART*UII - KFST4_PART
+   IF (STORE_RADIATION_TERMS) THEN
+      RADIATION_EMISSION   = KFST4_GAS
+      RADIATION_ABSORPTION = KAPPA_GAS*UII
+   ENDIF
 ENDIF
 
 ! Calculate the incoming radiative flux onto the solid particles

Source/read.f90

@@ -15966,6 +15966,7 @@ SUBROUTINE GET_QUANTITY_INDEX(SMOKEVIEW_LABEL,SMOKEVIEW_BAR_LABEL,OUTPUT_INDEX,O
                               Y_INDEX,Z_INDEX,PART_INDEX,DUCT_INDEX,NODE_INDEX,REAC_INDEX,MATL_INDEX,OUTTYPE, &
                               QUANTITY,QUANTITY2,SPEC_ID_IN,PART_ID,DUCT_ID,NODE_ID,REAC_ID,MATL_ID,CELL_L,&
                               DUCT_CELL_INDEX,SLICETYPE)
+
 CHARACTER(*), INTENT(INOUT) :: QUANTITY
 CHARACTER(*), INTENT(OUT) :: SMOKEVIEW_LABEL,SMOKEVIEW_BAR_LABEL
 CHARACTER(*) :: SPEC_ID_IN,PART_ID,DUCT_ID,NODE_ID
@@ -16165,6 +16166,8 @@ SUBROUTINE GET_QUANTITY_INDEX(SMOKEVIEW_LABEL,SMOKEVIEW_BAR_LABEL,OUTPUT_INDEX,O
     QUANTITY=='TOTAL MASS FLUX X'     .OR. QUANTITY=='TOTAL MASS FLUX Y'     .OR.  QUANTITY=='TOTAL MASS FLUX Z'     .OR. &
     QUANTITY=='TOTAL MASS FLUX WALL') STORE_SPECIES_FLUX = .TRUE.
 
+IF (QUANTITY=='RADIATION EMISSION' .OR. QUANTITY=='RADIATION ABSORPTION') STORE_RADIATION_TERMS = .TRUE.
+
 IF (TRIM(QUANTITY)=='HRRPUV REAC') THEN
    IF (TRIM(REAC_ID)=='null') THEN
       WRITE(MESSAGE,'(A)') 'ERROR(1011): HRRPUV REAC requires a REAC_ID.'