diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 7cf1a6ceed5..2e793180f71 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -10205,7 +10205,8 @@ \subsection{Heat Release Rate and Energy Conservation} \end{eqnarray} \begin{description} \item[{\ct Q\_ENTH}] The change in the sensible enthalpy of the gas. $\rho$ is the density of the gas (kg/m$^3$). $h_{\rm s}$ is the \underline{s}ensible enthalpy of the gas (kJ/kg). The volume integral is over the entire domain. -\item[{\ct HRR}] The total heat release rate of the fire (kW). By default the effects of any surface oxidation reactions are included in this value. This helps when comparing to heat release measurements obtained from oxygen consumption calorimetry, as the additional oxygen sink from surface reactions will be lumped into the measurement. However, it is possible to output only the contribution from gas-phase combustion by setting {\ct HRR\_GAS\_ONLY=T} on the {\ct DUMP} line. +\item[{\ct HRR}] The heat release rate of the fire (kW) resulting from gas phase combustion. +\item[{\ct HRR\_OX}] The heat release rate (kW) of any surface oxidation reactions. This helps when comparing to heat release measurements obtained from oxygen consumption calorimetry, as the additional oxygen sink from surface reactions will be lumped into the measurement. \item[{\ct Q\_RADI}] The thermal radiation {\em into} the domain from the exterior boundary or particles. $\dot{\bq}_{\rm r}''$ is the \underline{r}adiation heat flux vector (\unit{kW/m^2}). Its divergence represents the net radiative emission from a volume of gas. Typically, {\ct Q\_RADI} has a negative value, meaning that a fire or hot gases radiate energy out of the domain. $\dq_{\rm p,r}$ is the \underline{r}adiation absorbed by a droplet or \underline{p}article (kW). This term is added to {\ct Q\_RADI} and subtracted from {\ct Q\_PART} because it is implicitly included in $\nabla \cdot \dot{\bq}_{\rm r}''$ and needs to be separated off for the purpose of explicitly accounting for it in the energy budget. \item[{\ct Q\_CONV}] The flow of sensible enthalpy {\em into} the computational domain. $\dm_{\rm p,\alpha}$ is the production rate of gas species $\alpha$ from a solid \underline{p}article or liquid droplet (kg/s). $h_{\rm s,\alpha}$ is the \underline{s}ensible enthalpy of gas species $\alpha$ (kJ/kg). $\rho$ is the gas density (kg/m$^3$), $\bu$ is the velocity vector (m/s). $h_{\rm s}$ is the \underline{s}ensible enthalpy of the gas. If the gas is flowing out of the domain, $\bu \cdot \d {\bf S}$ is positive. \item[{\ct Q\_COND}] The convective heat flux {\em into} the computational domain. $\dq_{\rm c}''$ is the heat \underline{c}onvected from the gas to a surface. If the gas is relatively hot and the surfaces/particles/droplets relatively cool, {\ct Q\_COND} is negative. At {\ct OPEN} boundaries, {\ct Q\_COND} is $\int k \nabla T \cdot \d {\bf S}$, where $k$ (kW/(m$\cdot$K)) is the turbulent thermal conductivity of the gas and $\nabla T$ is the temperature gradient across the open boundary. $\dq_{\rm p,w}$ is the energy transferred from a solid surface (\underline{w}all) to a droplet or \underline{p}article adhering to it. Notice that it is subtracted off in {\ct Q\_PART} because it makes no contribution to the energy of the gas. @@ -12066,7 +12067,6 @@ \section{\texorpdfstring{{\tt DUMP}}{DUMP} (Output Parameters)} {\ct DT\_TMP} & Real & Section~\ref{info:CSVF} & s & \\ \hline {\ct DT\_UVW} & Real & Section~\ref{info:CSVF} & s & \\ \hline {\ct FLUSH\_FILE\_BUFFERS} & Logical & Section~\ref{info:DUMP} & & {\ct T} \\ \hline -{\ct HRR\_GAS\_ONLY} & Logical & Section~\ref{info:HRR} & & {\ct F} \\ \hline {\ct MASS\_FILE} & Logical & Section~\ref{info:DUMP} & & {\ct F} \\ \hline {\ct MAXIMUM\_PARTICLES} & Integer & Section~\ref{info:controlling_droplets}& & 1000000 \\ \hline {\ct NFRAMES} & Integer & Section~\ref{info:DUMP} & & 1000 \\ \hline diff --git a/Source/cons.f90 b/Source/cons.f90 index 2219180275c..c6d6506199c 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -260,7 +260,6 @@ MODULE GLOBAL_CONSTANTS LOGICAL :: PERIODIC_DOMAIN_Y=.FALSE. !< The domain is periodic \f$ y \f$ LOGICAL :: PERIODIC_DOMAIN_Z=.FALSE. !< The domain is periodic \f$ z \f$ LOGICAL :: OPEN_WIND_BOUNDARY=.FALSE. !< There is a prevailing wind -LOGICAL :: HRR_GAS_ONLY=.FALSE. !< Surface oxidation is not included in total HRR LOGICAL :: WRITE_DEVC_CTRL=.FALSE. !< Flag for writing DEVC and CTRL logfile LOGICAL :: INIT_INVOKED_BY_SURF=.FALSE. !< Flag indicating that a SURF line specifies an INIT line LOGICAL :: NO_PRESSURE_ZONES=.FALSE. !< Flag to suppress pressure zones diff --git a/Source/data.f90 b/Source/data.f90 index 628108ad649..33c8fdfffbe 100644 --- a/Source/data.f90 +++ b/Source/data.f90 @@ -7,7 +7,7 @@ MODULE OUTPUT_DATA IMPLICIT NONE (TYPE,EXTERNAL) -INTEGER, PARAMETER :: N_Q_DOT=8 +INTEGER, PARAMETER :: N_Q_DOT=9 INTEGER :: PLOT3D_QUANTITY_INDEX(5),PLOT3D_Y_INDEX(5)=0,PLOT3D_Z_INDEX(5)=0,PLOT3D_PART_INDEX(5),& PLOT3D_VELO_INDEX(5)=0 CHARACTER(LABEL_LENGTH) :: PLOT3D_QUANTITY(5),PLOT3D_SPEC_ID(5),PLOT3D_PART_ID(5),PLOT3D_SMOKEVIEW_BAR_LABEL(5) diff --git a/Source/dump.f90 b/Source/dump.f90 index f28354f3c7c..1fe69412445 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -799,13 +799,13 @@ SUBROUTINE INITIALIZE_GLOBAL_DUMPS(T,DT) CALL APPEND_FILE(LU_HRR,2,T_BEGIN+(T-T_BEGIN)*TIME_SHRINK_FACTOR) ELSE OPEN(LU_HRR,FILE=FN_HRR,FORM='FORMATTED',STATUS='REPLACE') - WRITE(TCFORM,'(A,I0,A)') "(",9+N_TRACKED_SPECIES+N_ZONE_TMP,"(A,','),A)" - WRITE(LU_HRR,TCFORM) 's','kW','kW','kW','kW','kW','kW','kW','kW','kW',('kg/s',N=1,N_TRACKED_SPECIES),('Pa',N=1,N_ZONE_TMP) + WRITE(TCFORM,'(A,I0,A)') "(",N_Q_DOT+1+N_TRACKED_SPECIES+N_ZONE_TMP,"(A,','),A)" + WRITE(LU_HRR,TCFORM) 's','kW','kW','kW','kW','kW','kW','kW','kW','kW','kW',('kg/s',N=1,N_TRACKED_SPECIES),('Pa',N=1,N_ZONE_TMP) IF (N_ZONE_TMP>0) THEN - WRITE(LU_HRR,TCFORM) 'Time','HRR','Q_RADI','Q_CONV','Q_COND','Q_DIFF','Q_PRES','Q_PART','Q_ENTH','Q_TOTAL',& + WRITE(LU_HRR,TCFORM) 'Time','HRR','HRR_OX','Q_RADI','Q_CONV','Q_COND','Q_DIFF','Q_PRES','Q_PART','Q_ENTH','Q_TOTAL',& ('MLR_'//TRIM(SPECIES_MIXTURE(N)%ID),N=1,N_TRACKED_SPECIES),(TRIM(P_ZONE(N)%ID),N=1,N_ZONE_TMP) ELSE - WRITE(LU_HRR,TCFORM) 'Time','HRR','Q_RADI','Q_CONV','Q_COND','Q_DIFF','Q_PRES','Q_PART','Q_ENTH','Q_TOTAL',& + WRITE(LU_HRR,TCFORM) 'Time','HRR','HRR_OX','Q_RADI','Q_CONV','Q_COND','Q_DIFF','Q_PRES','Q_PART','Q_ENTH','Q_TOTAL',& ('MLR_'//TRIM(SPECIES_MIXTURE(N)%ID),N=1,N_TRACKED_SPECIES) ENDIF ENDIF @@ -3923,7 +3923,7 @@ SUBROUTINE WRITE_DIAGNOSTICS(T,DT) IF (SIM_MODE==DNS_MODE .OR. CHECK_VN) WRITE(LU_OUTPUT,230) M%VN,M%I_VN,M%J_VN,M%K_VN IF (M%NLP>0) WRITE(LU_OUTPUT,141) M%NLP IF (ABS(Q_DOT(1,NM))>1._EB) WRITE(LU_OUTPUT,119) Q_DOT(1,NM)/1000._EB - IF (ABS(Q_DOT(2,NM))>1._EB) WRITE(LU_OUTPUT,120) Q_DOT(2,NM)/1000._EB + IF (ABS(Q_DOT(3,NM))>1._EB) WRITE(LU_OUTPUT,120) Q_DOT(3,NM)/1000._EB IF (M%DT_RESTRICT_STORE>0 ) THEN WRITE(LU_OUTPUT,121) M%DT_RESTRICT_STORE M%DT_RESTRICT_STORE=0 @@ -9935,13 +9935,14 @@ END SUBROUTINE DUMP_HVAC !> \param NM Mesh number !> \details !> Q_DOT(1,NM) = \f$ \int \dot{q}''' \, dV \f$ -!> Q_DOT(2,NM) = \f$ \int \nabla \cdot \mathbf{q}_{\rm r}'' \, dV \f$ -!> Q_DOT(3,NM) = \f$ \int \mathbf{u} \rho h_{\rm s} \cdot \, d\mathbf{S} \f$ -!> Q_DOT(4,NM) = \f$ \int k \nabla T \cdot d\mathbf{S} \f$ -!> Q_DOT(5,NM) = \f$ \int \sum_\alpha h_{{\rm s},\alpha} \rho D_\alpha \nabla Z_\alpha \cdot d\mathbf{S} \f$ -!> Q_DOT(6,NM) = \f$ \int dp/dt \, dV \f$ -!> Q_DOT(7,NM) = \f$ \sum \dot{q}_{\rm p} \f$ -!> Q_DOT(8,NM) = \f$ \int d(\rho h_{\rm s})/dt \, dV \f$ +!> Q_DOT(2,NM) = \f$ \int \dot{q}_{\rm ox}'' \, dS \f$ +!> Q_DOT(3,NM) = \f$ \int \nabla \cdot \mathbf{q}_{\rm r}'' \, dV \f$ +!> Q_DOT(4,NM) = \f$ \int \mathbf{u} \rho h_{\rm s} \cdot \, d\mathbf{S} \f$ +!> Q_DOT(5,NM) = \f$ \int k \nabla T \cdot d\mathbf{S} \f$ +!> Q_DOT(6,NM) = \f$ \int \sum_\alpha h_{{\rm s},\alpha} \rho D_\alpha \nabla Z_\alpha \cdot d\mathbf{S} \f$ +!> Q_DOT(7,NM) = \f$ \int dp/dt \, dV \f$ +!> Q_DOT(8,NM) = \f$ \sum \dot{q}_{\rm p} \f$ +!> Q_DOT(9,NM) = \f$ \int d(\rho h_{\rm s})/dt \, dV \f$ SUBROUTINE UPDATE_HRR(DT,NM) @@ -9973,8 +9974,8 @@ SUBROUTINE UPDATE_HRR(DT,NM) IF (CYLINDRICAL) VC = VC*2._EB*PI Q_DOT(1,NM) = Q_DOT(1,NM) + Q(I,J,K)*VC - Q_DOT(2,NM) = Q_DOT(2,NM) + QR(I,J,K)*VC - Q_DOT(6,NM) = Q_DOT(6,NM) + 0.5_EB*(D_PBAR_DT_S(PRESSURE_ZONE(I,J,K))+D_PBAR_DT(PRESSURE_ZONE(I,J,K)))*VC + Q_DOT(3,NM) = Q_DOT(3,NM) + QR(I,J,K)*VC + Q_DOT(7,NM) = Q_DOT(7,NM) + 0.5_EB*(D_PBAR_DT_S(PRESSURE_ZONE(I,J,K))+D_PBAR_DT(PRESSURE_ZONE(I,J,K)))*VC ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES) CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMP(I,J,K)) ENTHALPY_SUM(NM) = ENTHALPY_SUM(NM) + RHO(I,J,K)*H_S*VC @@ -9982,12 +9983,12 @@ SUBROUTINE UPDATE_HRR(DT,NM) ENDDO ENDDO -IF (CC_IBM) CALL ADD_Q_DOT_CUTCELLS(NM,Q_DOT(1,NM),Q_DOT(2,NM),Q_DOT(6,NM),ENTHALPY_SUM(NM)) +IF (CC_IBM) CALL ADD_Q_DOT_CUTCELLS(NM,Q_DOT(1,NM),Q_DOT(3,NM),Q_DOT(7,NM),ENTHALPY_SUM(NM)) IF (ICYC>0) THEN - Q_DOT(8,NM) = (ENTHALPY_SUM(NM)-ENTHALPY_SUM_OLD)/DT + Q_DOT(9,NM) = (ENTHALPY_SUM(NM)-ENTHALPY_SUM_OLD)/DT ELSE - Q_DOT(8,NM) = 0._EB + Q_DOT(9,NM) = 0._EB ENDIF ! Compute the surface integral of all Del Dot terms @@ -10036,10 +10037,10 @@ SUBROUTINE UPDATE_HRR(DT,NM) AREA_F = B1%AREA IF (TWO_D) AREA_F = AREA_F/DY(BC%JJG) IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI - Q_DOT(3,NM) = Q_DOT(3,NM) - U_N*B1%RHO_F*H_S*AREA_F - Q_DOT(4,NM) = Q_DOT(4,NM) - B1%Q_CON_F*AREA_F - Q_DOT(5,NM) = Q_DOT(5,NM) - H_S_J_ALPHA*AREA_F - IF (.NOT. HRR_GAS_ONLY) Q_DOT(1,NM) = Q_DOT(1,NM) + B1%Q_DOT_O2_PP*AREA_F + Q_DOT(2,NM) = Q_DOT(2,NM) + B1%Q_DOT_O2_PP*AREA_F + Q_DOT(4,NM) = Q_DOT(4,NM) - U_N*B1%RHO_F*H_S*AREA_F + Q_DOT(5,NM) = Q_DOT(5,NM) - B1%Q_CON_F*AREA_F + Q_DOT(6,NM) = Q_DOT(6,NM) - H_S_J_ALPHA*AREA_F ENDDO WALL_LOOP CFACE_LOOP : DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS @@ -10068,13 +10069,13 @@ SUBROUTINE UPDATE_HRR(DT,NM) AREA_F = B1%AREA IF (TWO_D) AREA_F = AREA_F/DY(BC%JJG) IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI - Q_DOT(3,NM) = Q_DOT(3,NM) - U_N*B1%RHO_F*H_S*AREA_F - Q_DOT(4,NM) = Q_DOT(4,NM) - B1%Q_CON_F*AREA_F - Q_DOT(5,NM) = Q_DOT(5,NM) - H_S_J_ALPHA*AREA_F + Q_DOT(4,NM) = Q_DOT(4,NM) - U_N*B1%RHO_F*H_S*AREA_F + Q_DOT(5,NM) = Q_DOT(5,NM) - B1%Q_CON_F*AREA_F + Q_DOT(6,NM) = Q_DOT(6,NM) - H_S_J_ALPHA*AREA_F ENDDO CFACE_LOOP -IF (OXIDATION_REACTION .AND. .NOT. HRR_GAS_ONLY) THEN +IF (OXIDATION_REACTION) THEN PARTICLE_LOOP: DO IP=1,NLP LP => LAGRANGIAN_PARTICLE(IP) LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX) @@ -10084,7 +10085,7 @@ SUBROUTINE UPDATE_HRR(DT,NM) AREA_F = B1%AREA IF (TWO_D) AREA_F = AREA_F/DY(BC%JJG) IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI - Q_DOT(1,NM) = Q_DOT(1,NM) + LP%PWT*B1%Q_DOT_O2_PP*AREA_F + Q_DOT(2,NM) = Q_DOT(2,NM) + LP%PWT*B1%Q_DOT_O2_PP*AREA_F ENDDO PARTICLE_LOOP ENDIF @@ -10152,7 +10153,7 @@ SUBROUTINE DUMP_HRR(T,DT) ENDDO ENDIF -WRITE(TCFORM,'(A,I0,5A)') "(",9+N_TRACKED_SPECIES+N_ZONE_TMP,"(",FMT_R,",','),",FMT_R,")" +WRITE(TCFORM,'(A,I0,5A)') "(",N_Q_DOT+1+N_TRACKED_SPECIES+N_ZONE_TMP,"(",FMT_R,",','),",FMT_R,")" IF (N_ZONE_TMP>0) THEN WRITE(LU_HRR,TCFORM) STIME,0.001_EB*Q_DOT_TOTAL(1:N_Q_DOT),0.001_EB*SUM(Q_DOT_TOTAL(1:N_Q_DOT-1)),& M_DOT_TOTAL(1:N_TRACKED_SPECIES),(P_ZONE_P(I),I=1,N_ZONE_TMP) diff --git a/Source/part.f90 b/Source/part.f90 index 4ed9e0465bf..e76bd3e4540 100644 --- a/Source/part.f90 +++ b/Source/part.f90 @@ -3907,9 +3907,9 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) ! Add energy losses and gains to overall energy budget array - Q_DOT(7,NM) = Q_DOT(7,NM) - Q_RAD*WGT/DT ! Q_PART - Q_DOT(3,NM) = Q_DOT(3,NM) + M_VAP*H_S_B*WGT/DT ! Q_CONV - Q_DOT(2,NM) = Q_DOT(2,NM) + Q_RAD*WGT/DT ! Q_RADI + Q_DOT(8,NM) = Q_DOT(8,NM) - Q_RAD*WGT/DT ! Q_PART + Q_DOT(4,NM) = Q_DOT(4,NM) + M_VAP*H_S_B*WGT/DT ! Q_CONV + Q_DOT(3,NM) = Q_DOT(3,NM) + Q_RAD*WGT/DT ! Q_RADI IF (LPC%Z_INDEX>0) M_DOT(LPC%Z_INDEX,NM) = M_DOT(LPC%Z_INDEX,NM) + WGT*M_VAP/DT/LPC%ADJUST_EVAPORATION @@ -4282,10 +4282,10 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) ! Add energy losses and gains to overall energy budget array - Q_DOT(7,NM) = Q_DOT(7,NM) - (Q_CON_GAS + Q_CON_WALL + Q_RAD)*WGT/DT ! Q_PART - Q_DOT(3,NM) = Q_DOT(3,NM) + M_VAP*H_S_B*WGT/DT ! Q_CONV - Q_DOT(2,NM) = Q_DOT(2,NM) + Q_RAD*WGT/DT ! Q_RADI - Q_DOT(4,NM) = Q_DOT(4,NM) + Q_CON_WALL*WGT/DT ! Q_COND + Q_DOT(8,NM) = Q_DOT(8,NM) - (Q_CON_GAS + Q_CON_WALL + Q_RAD)*WGT/DT ! Q_PART + Q_DOT(4,NM) = Q_DOT(4,NM) + M_VAP*H_S_B*WGT/DT ! Q_CONV + Q_DOT(3,NM) = Q_DOT(3,NM) + Q_RAD*WGT/DT ! Q_RADI + Q_DOT(5,NM) = Q_DOT(5,NM) + Q_CON_WALL*WGT/DT ! Q_COND IF (LPC%Z_INDEX>0) M_DOT(LPC%Z_INDEX,NM) = M_DOT(LPC%Z_INDEX,NM) + WGT*M_VAP/DT/LPC%ADJUST_EVAPORATION diff --git a/Source/read.f90 b/Source/read.f90 index 13eb5f04a80..f6b764a1706 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -2322,7 +2322,7 @@ SUBROUTINE READ_DUMP DIAGNOSTICS_INTERVAL,& DT_BNDF,DT_CPU,DT_CTRL,DT_DEVC,DT_FLUSH,DT_HRR,DT_HVAC,DT_ISOF,DT_MASS,DT_PART,DT_PL3D,DT_PROF,& DT_RADF,DT_RESTART,DT_SL3D,DT_SLCF,DT_SMOKE3D,DT_UVW,DT_TMP,DT_SPEC,& - FLUSH_FILE_BUFFERS,GET_CUTCELLS_VERBOSE,HRR_GAS_ONLY,MASS_FILE,MAXIMUM_PARTICLES,MMS_TIMER,& + FLUSH_FILE_BUFFERS,GET_CUTCELLS_VERBOSE,MASS_FILE,MAXIMUM_PARTICLES,MMS_TIMER,& NFRAMES,PLOT3D_PART_ID,PLOT3D_QUANTITY,PLOT3D_SPEC_ID,PLOT3D_VELO_INDEX,& RAMP_BNDF,RAMP_CPU,RAMP_CTRL,RAMP_DEVC,RAMP_FLUSH,RAMP_HRR,RAMP_HVAC,RAMP_ISOF,RAMP_MASS,& RAMP_PART,RAMP_PL3D,RAMP_PROF,RAMP_RADF,RAMP_RESTART,RAMP_SLCF,RAMP_SL3D,RAMP_SMOKE3D,& diff --git a/Source/wall.f90 b/Source/wall.f90 index 0c10fa1ed16..bee06939fc5 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -1425,7 +1425,7 @@ SUBROUTINE DEPOSIT_PARTICLE_MASS(NM,LP,LPC) ZZ_GET(NS) = 1._EB CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S_B,B1%TMP_F) !$OMP CRITICAL - Q_DOT(3,NM) = Q_DOT(3,NM) + B1%M_DOT_G_PP_ADJUST(NS)*B1%AREA*H_S_B*LP%PWT ! Q_CONV + Q_DOT(4,NM) = Q_DOT(4,NM) + B1%M_DOT_G_PP_ADJUST(NS)*B1%AREA*H_S_B*LP%PWT ! Q_CONV !$OMP END CRITICAL ENDDO @@ -1444,8 +1444,8 @@ SUBROUTINE DEPOSIT_PARTICLE_MASS(NM,LP,LPC) ! Add energy losses and gains to overall energy budget array -Q_DOT(7,NM) = Q_DOT(7,NM) - (B1%Q_CON_F + B1%Q_RAD_IN - B1%Q_RAD_OUT)*B1%AREA*LP%PWT ! Q_PART -Q_DOT(2,NM) = Q_DOT(2,NM) + (B1%Q_RAD_IN-B1%Q_RAD_OUT)*B1%AREA*LP%PWT ! Q_RADI +Q_DOT(8,NM) = Q_DOT(8,NM) - (B1%Q_CON_F + B1%Q_RAD_IN - B1%Q_RAD_OUT)*B1%AREA*LP%PWT ! Q_PART +Q_DOT(3,NM) = Q_DOT(3,NM) + (B1%Q_RAD_IN-B1%Q_RAD_OUT)*B1%AREA*LP%PWT ! Q_RADI !$OMP END CRITICAL ! Calculate the mass flux of fuel gas from particles