FireX: Merge with firemodels/master

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

@@ -5025,24 +5025,99 @@ \section{SFPE Verification Cases}
 \subsection{Case 6: 2-D Heat Transfer with Cooling by Convection}
 \label{SFPE_Case_6}
 
-A 2~m by 2~m square column ($k=1$~W/(m~K), $\rho=1$~kg/m$^3$, $c=0.001$~kJ/(kg~K), $\epsilon=0$) with an initial temperature of 1000~°C cools via convection only. Assuming that $h=1$~W/(m$^2$~K) and the surrounding air temperature is 0~°C, calculate the temperature at the center of the column as a function of time.
+A 2~m by 2~m square column ($k=1$~W/(m~K), $\rho=1$~kg/m$^3$, $c=0.001$~kJ/(kg~K), $\epsilon=0$) with an initial temperature of 1000~°C cools via convection only. Assuming that $h=1$~W/(m$^2$~K) and the surrounding air temperature is 0~°C, calculate the temperature at the center of the column as a function of time (Fig.~\ref{fig:SFPE_Case_6}).
 
-This is an unusual set of parameters, but nevertheless this case is used to test a multi-dimensional heat transfer solver. The results are shown in Fig.~\ref{fig:SFPE_Case_6_7}.
+\begin{figure}[ht]
+\centering
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_6}
+\caption[The SFPE heat transfer verification Case 6]{Temperature at the center of a 2~m by 2~m column that is initially hot and is cooled via convection by surrounding air.}
+\label{fig:SFPE_Case_6}
+\end{figure}
+
+\FloatBarrier
 
 \subsection{Case 7: 2-D Heat Transfer by Convection and Radiation}
 \label{SFPE_Case_7}
 
 A 0.2 m by 0.2 m square column ($k=1$~W/(m~K), $\rho=2400$~kg/m$^3$, $c=1$~kJ/(kg~K), $\epsilon=0.8$) is heated according to the ISO~834 time-temperature curve
 \be
-   T(t) = T_\infty + 345 \, \ln \left( 8t/60 + 1 \right)
+   T(t) = T_\infty + 345 \, \ln \left( 8t/60 + 1 \right) \label{ISO_834}
 \ee
-where the time, $t$, is in seconds. Assuming that $h=10$~W/(m$^2$~K) and that the initial temperature is $T_\infty=273$~K, calculate the temperature at the column center, corner and middle side surface as a function of time. The results are shown in Fig.~\ref{fig:SFPE_Case_6_7}.
+where the time, $t$, is in seconds. Assuming that $h=10$~W/(m$^2$~K) and that the initial temperature is $T_\infty=273$~K, calculate the temperature at the column center, corner and middle side surface as a function of time (Fig.~\ref{fig:SFPE_Case_7}).
 
 \begin{figure}[ht]
-\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_6}
+\centering
 \includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_7}
-\caption[The SFPE heat transfer verification cases 6 and 7]{Comparisons of FDS 2-D heat transfer calculations with solutions generated by a finite-element solver.}
-\label{fig:SFPE_Case_6_7}
+\caption[The SFPE heat transfer verification Case 7]{Temperature at the center, corner and side of a 20~cm by 20~cm column heated by a fire that follows a specified time-temperature curve.}
+\label{fig:SFPE_Case_7}
+\end{figure}
+
+\FloatBarrier
+
+
+\subsection{Case 8: 2-D Heat Transfer with Temperature-Dependent Conductivity}
+\label{SFPE_Case_8}
+
+A 0.2~m by 0.2~m square column ($\rho=2400$~kg/m$^3$, $c=1$~kJ/(kg~K), $\epsilon=0.8$) is heated according to the ISO~834 time-temperature curve, Eq.~(\ref{ISO_834}. The thermal conductivity of the column material varies linearly with temperature such that its value is 1.5~W/(m~K) at 0~°C, 0.7~W/(m~K) at 200~°C, and 0.5~W/(m~K) at 1000~°C. Assuming that $h=10$~W/(m$^2$~K) and that the initial air temperature is 0~°C, calculate the temperature at the column center, corner and middle side surface as a function of time (Fig.~\ref{fig:SFPE_Case_8}).
+
+\begin{figure}[ht]
+\centering
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_8}
+\caption[The SFPE heat transfer verification Case 8]{Temperature at the center, corner and side of a 20~cm by 20~cm column heated by a fire that follows a specified time-temperature curve. The thermal conductivity of the material varies with temperature.}
+\label{fig:SFPE_Case_8}
+\end{figure}
+
+\FloatBarrier
+
+
+\subsection{Case 9: 2-D Heat Transfer in a Composite Section with Temperature-Dependent Conductivity}
+\label{SFPE_Case_9}
+
+A hollow square metal tube ($\rho=7850$~kg/m$^3$, $c=0.6$~kJ/(kg~K), $\epsilon=0.8$) is filled with an insulation material ($k=0.05$~W/(m~K), $\rho=50$~kg/m$^3$, $c=1$~kJ/(kg~K)). The thermal conductivity of the metal tube varies linearly with temperature such that its value is 54.7~W/(m~K) at 0~°C, 27.3~W/(m~K) at 800~°C, and 27.3~W/(m~K) at 1200~°C. The tube walls are 0.5~mm thick, and the exterior width of the assembly is 0.201~m. The surrounding air temperature is 1000~°C, and the initial temperature of the assembly is 0~°C. Assuming that the heating is by convection and radiation, and that the heat transfer coefficient is 10~W/(m$^2$~K), calculate the temperature at the center of the tube as a function of time (Fig.~\ref{fig:SFPE_Case_9}).
+
+Normally, a 3-D conducting solid with a thin ``skin'' would be modeled by assigning a thin layer to the \ct{SURF} line that describes the surface properties of the obstruction. However, in this case, the difference in the enthalpy, $c_{\rm s} \rho_{\rm s}$, between the metal and the insulation material leads to inaccuracies in the alternating direction implicit (ADI) solver. In short, the energy absorbed by the thin metal skin is transferred from the $x$ or $y$ direction to the $z$ direction, which has no metal skin. Thus, the temperature in the $z$ direction becomes fictitiously high because its thermal mass is relatively low. As a work-around, the metal skin is thickened from 0.5~mm to 10~mm, the width of a gas phase cell, while its density is decreased by a factor of 20.
+
+\begin{figure}[ht]
+\centering
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_9}
+\caption[The SFPE heat transfer verification Case 9]{Temperature at the center of a 20~cm by 20~cm thin metal tube filled with insulation material.}
+\label{fig:SFPE_Case_9}
+\end{figure}
+
+\FloatBarrier
+
+\subsection{Case 13: 2-D Heat Transfer with Moisture Evaporation}
+\label{SFPE_Case_13}
+
+A 0.2~m by 0.2~m square column (dry properties: $\rho_{\rm d}=2400$~kg/m$^3$, $c_{\rm d}=1$~kJ/(kg~K), $\epsilon=0.8$) is heated according to the ISO~834 time-temperature curve, Eq.~(\ref{ISO_834}). The thermal conductivity of the column material varies linearly with temperature such that its value is 1.5~W/(m~K) at 0~°C, 0.7~W/(m~K) at 200~°C, and 0.5~W/(m~K) at 1000~°C. The mass fraction of water within the column is $Y_{\rm w}=0.0208$, and it is assumed that the water evaporates at temperatures between 100~°C and 120~°C. The density of water and specific heat capacity of water can be taken as $\rho_{\rm w}=1000$~kg/m$^3$ and $c_{\rm w}=4.187$~kJ/(kg~K), respectively. The latent heat of evaporation of water ($h_{\rm v}=2260$~kJ/kg) is assumed to be added to the specific heat of the material. During evaporation, the amount of water is assumed to decrease linearly to zero. Assuming that $h=10$~W/(m$^2$~K) and that the initial temperature is 20~°C, calculate the temperature at the column center, corner and middle side surface as a function of time (Fig.~\ref{fig:SFPE_Case_13}).
+
+Normally, this scenario would be modeled in FDS by assuming that the water boils in the neighborhood of 100~$^\circ$C following an Arhennius reaction. However, this method would not allow for a linear decrease in the amount of water during its evaporation from 100~$^\circ$C to 120~$^\circ$C. Instead, the latent heat of vaporization of water has been incorporated into an effective specific heat of the column material, shown in the left plot of Fig.~\ref{fig:SFPE_Case_13}. The area of the triangular region is the product of the latent heat of evaporation and the mass fraction of water, $h_{\rm v} Y_{\rm w}$. The effective specific heat of the wet column is given by
+\be
+   c_{\rm eff} = \frac{ m_{\rm d}''' c_{\rm d} + m_{\rm w}''' c_{\rm w} }{ \rho_{\rm d} } \approx 1.089 \; \hbox{kJ/(kg K)} \quad ; \quad m_{\rm d}''' = 2400 \; \hbox{kg/m}^3 \quad ; \quad m_{\rm w}''' = \frac{ Y_{\rm w}}{1-Y_{\rm w}} \,  m_{\rm d}''' \approx 50.1  \hbox{kg/m}^3
+\ee
+where $m_{\rm d}'''$ is the mass per unit volume of the dry column material and $m_{\rm w}'''$ is the mass per unit volume of the water that fills the pores of the column material. Note that it is assumed that the volume of the wet column is the same as that of the dry. Also note that FDS does not allow for a change in material density without a reaction. Thus, in this case, the density of the column is taken as its dry density, which is seen in the denominator of the formula for the effective specific heat of the wet column.
+
+\begin{figure}[ht]
+\centering
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_13b}
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_13}
+\caption[The SFPE heat transfer verification Case 13]{(Left) The effective specific heat of the column material that incorporates the latent heat of evaporation of water. (Right) The temperature at the center, corner and side of a 20~cm by 20~cm column heated by a fire that follows a specified time-temperature curve.}
+\label{fig:SFPE_Case_13}
+\end{figure}
+
+
+\FloatBarrier
+
+\subsection{Case 16: 3-D Heat Transfer with Non-Uniform Heat Flux}
+\label{SFPE_Case_16}
+
+A 1~m by 2~m rectangular metal plate ($k=50$~W/(m~K), $\rho=7850$~kg/m$^3$, $c=0.5$~kJ/(kg~K), $\epsilon=0.8$) that is 10~cm thick is heated over a quarter of its upper surface with an incident radiant heat flux of 30~kW/m$^2$ with convection to ambient. In Cartesian coordinates $(x,y)$, the corners of the plate are $(0,0)$, $(1,0)$, $(1,2)$, and $(0,2)$. All lengths are in meters. The heating occurs at the top of the plate within the rectangular area whose corners are $(0,0)$, $(0.5,0)$, $(0.5,1)$, $(0,1)$. The bottom surface and the remainder of the top surface are cooled by convection and radiation to ambient. The sides of the plate are insulated. Assuming that the surrounding air temperature is 20~°C, and that the heat transfer coefficient is 10~W/(m$^2$~K), calculate the temperature along the lines $x=0.25$~m and $y=0.5$~m at the mid-depth of the plate at 60~min (Fig.~\ref{fig:SFPE_Case_16}).
+
+\begin{figure}[ht]
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_16a}
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_16b}
+\caption[The SFPE heat transfer verification Case 16]{Temperature along the length (left) and width (right) of a metal plate heated over a quarter of its surface by a radiant heater of 30~kW/m$^2$ for 60~min.}
+\label{fig:SFPE_Case_16}
 \end{figure}
 
 

Source/data.f90

@@ -1280,21 +1280,21 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(534)%SHORT_NAME = 'rho*(Z*u + D*dZ/dx)'
 OUTPUT_QUANTITY(534)%CELL_POSITION = CELL_FACE
 OUTPUT_QUANTITY(534)%IOR = 1
-OUTPUT_QUANTITY(534)%SPEC_ID_REQUIRED=.TRUE.
+!OUTPUT_QUANTITY(534)%SPEC_ID_REQUIRED=.TRUE.
 
 OUTPUT_QUANTITY(535)%NAME = 'TOTAL MASS FLUX Y'
 OUTPUT_QUANTITY(535)%UNITS = 'kg/s/m2'
 OUTPUT_QUANTITY(535)%SHORT_NAME = 'rho*(Z*v + D*dZ/dy)'
 OUTPUT_QUANTITY(535)%CELL_POSITION = CELL_FACE
 OUTPUT_QUANTITY(535)%IOR = 2
-OUTPUT_QUANTITY(535)%SPEC_ID_REQUIRED=.TRUE.
+!OUTPUT_QUANTITY(535)%SPEC_ID_REQUIRED=.TRUE.
 
 OUTPUT_QUANTITY(536)%NAME = 'TOTAL MASS FLUX Z'
 OUTPUT_QUANTITY(536)%UNITS = 'kg/s/m2'
 OUTPUT_QUANTITY(536)%SHORT_NAME = 'rho*(Z*w + D*dZ/dz)'
 OUTPUT_QUANTITY(536)%CELL_POSITION = CELL_FACE
 OUTPUT_QUANTITY(536)%IOR = 3
-OUTPUT_QUANTITY(536)%SPEC_ID_REQUIRED=.TRUE.
+!OUTPUT_QUANTITY(536)%SPEC_ID_REQUIRED=.TRUE.
 
 OUTPUT_QUANTITY(550)%NAME = 'CUTCELL VELOCITY DIVERGENCE'
 OUTPUT_QUANTITY(550)%UNITS = '1/s'

Source/dump.f90

@@ -9657,29 +9657,40 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
       ENDIF
    CASE(534)  ! TOTAL MASS FLUX X
       GAS_PHASE_OUTPUT_RES = 0._EB
-      IF (Z_INDEX>0) THEN
-         GAS_PHASE_OUTPUT_RES = ADV_FX(II,JJ,KK,Z_INDEX) + DIF_FX(II,JJ,KK,Z_INDEX)
-      ELSEIF (Y_INDEX>0) THEN
-         GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
-                                (ADV_FX(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FX(II,JJ,KK,1:N_TRACKED_SPECIES)) )
+      IF (Z_INDEX<=0 .AND. Y_INDEX<=0) THEN
+         GAS_PHASE_OUTPUT_RES = SUM(ADV_FX(II,JJ,KK,:) + DIF_FX(II,JJ,KK,:))
+      ELSE
+         IF (Z_INDEX>0) THEN
+            GAS_PHASE_OUTPUT_RES = ADV_FX(II,JJ,KK,Z_INDEX) + DIF_FX(II,JJ,KK,Z_INDEX)
+         ELSEIF (Y_INDEX>0) THEN
+            GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
+                                   (ADV_FX(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FX(II,JJ,KK,1:N_TRACKED_SPECIES)) )
+         ENDIF
       ENDIF
    CASE(535)  ! TOTAL MASS FLUX Y
       GAS_PHASE_OUTPUT_RES = 0._EB
-      IF (Z_INDEX>0) THEN
-         GAS_PHASE_OUTPUT_RES = ADV_FY(II,JJ,KK,Z_INDEX) + DIF_FY(II,JJ,KK,Z_INDEX)
-      ELSEIF (Y_INDEX>0) THEN
-         GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
-                                (ADV_FY(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FY(II,JJ,KK,1:N_TRACKED_SPECIES)) )
+      IF (Z_INDEX<=0 .AND. Y_INDEX<=0) THEN
+         GAS_PHASE_OUTPUT_RES = SUM(ADV_FY(II,JJ,KK,:) + DIF_FY(II,JJ,KK,:))
+      ELSE
+         IF (Z_INDEX>0) THEN
+            GAS_PHASE_OUTPUT_RES = ADV_FY(II,JJ,KK,Z_INDEX) + DIF_FY(II,JJ,KK,Z_INDEX)
+         ELSEIF (Y_INDEX>0) THEN
+            GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
+                                   (ADV_FY(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FY(II,JJ,KK,1:N_TRACKED_SPECIES)) )
+         ENDIF
       ENDIF
    CASE(536)  ! TOTAL MASS FLUX Z
       GAS_PHASE_OUTPUT_RES = 0._EB
-      IF (Z_INDEX>0) THEN
-         GAS_PHASE_OUTPUT_RES = ADV_FZ(II,JJ,KK,Z_INDEX) + DIF_FZ(II,JJ,KK,Z_INDEX)
-      ELSEIF (Y_INDEX>0) THEN
-         GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
-                                (ADV_FZ(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FZ(II,JJ,KK,1:N_TRACKED_SPECIES)) )
+      IF (Z_INDEX<=0 .AND. Y_INDEX<=0) THEN
+         GAS_PHASE_OUTPUT_RES = SUM(ADV_FZ(II,JJ,KK,:) + DIF_FZ(II,JJ,KK,:))
+      ELSE
+         IF (Z_INDEX>0) THEN
+            GAS_PHASE_OUTPUT_RES = ADV_FZ(II,JJ,KK,Z_INDEX) + DIF_FZ(II,JJ,KK,Z_INDEX)
+         ELSEIF (Y_INDEX>0) THEN
+            GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
+                                   (ADV_FZ(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FZ(II,JJ,KK,1:N_TRACKED_SPECIES)) )
+         ENDIF
       ENDIF
-
     CASE(550) ! CUTCELL VELOCITY DIVERGENCE
        GAS_PHASE_OUTPUT_RES = CARTVELDIV(II,JJ,KK)
 

Source/pres.f90

@@ -2896,17 +2896,23 @@ SUBROUTINE ULMAT_H_MATRIX_SOLVER_SETUP(NM,IPZ)
    CALL HYPRE_IJMATRIXSETOBJECTTYPE(ZM%HYPRE_ZM%A_H, HYPRE_PARCSR, HYPRE_IERR)
    CALL HYPRE_IJMATRIXINITIALIZE_V2(ZM%HYPRE_ZM%A_H, HYPRE_MEMORY_HOST, HYPRE_IERR)
    IF (ZM%MTYPE==SYMM_INDEFINITE) THEN
-      ! Rows 1 to ZM%NUNKH-1, last column, set all to zero:
-      DO IROW=1,ZM%NUNKH-1
+      IF(ZM%NUNKH==1) THEN ! Single unknown, zero coefficient matrix (1,1). Set coefficient to 1.
+         NNZ_H_MAT(1) = 1
+         DEALLOCATE(JD_H_MAT); ALLOCATE(JD_H_MAT(1,1)); JD_H_MAT(1,1) = 1
+         DEALLOCATE( D_H_MAT); ALLOCATE( D_H_MAT(1,1));  D_H_MAT(1,1) = 1._EB
+      ELSE ! More than one unknown
+         ! Rows 1 to ZM%NUNKH-1, last column, set all to zero:
+         DO IROW=1,ZM%NUNKH-1
+            DO JCOL=1,NNZ_H_MAT(IROW)
+               IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) CYCLE ! Make zero matrix entries in last column.
+               D_H_MAT(JCOL,IROW) = 0._EB
+            ENDDO
+         ENDDO
+         ! Last row, all zeros except the diagonal that keeps diagonal number: Note after previous loop IROW==ZM%NUNKH
          DO JCOL=1,NNZ_H_MAT(IROW)
-            IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) CYCLE ! Make zero matrix entries in last column.
-            D_H_MAT(JCOL,IROW) = 0._EB
+            IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) D_H_MAT(JCOL,IROW) = 0._EB
          ENDDO
-      ENDDO
-      ! Last row, all zeros except the diagonal that keeps diagonal number: Note after previous loop IROW==ZM%NUNKH
-      DO JCOL=1,NNZ_H_MAT(IROW)
-         IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) D_H_MAT(JCOL,IROW) = 0._EB
-      ENDDO
+      ENDIF
    ENDIF
    JD_H_MAT = JD_H_MAT - 1
    DO IROW=1,ZM%NUNKH

Source/read.f90

@@ -14515,7 +14515,6 @@ SUBROUTINE READ_CTRL
          CF%CONTROL_INDEX = OR_GATE
       CASE('ONLY')
          CF%CONTROL_INDEX = XOR_GATE
-         CF%N = 1
       CASE('AT_LEAST')
          CF%CONTROL_INDEX = X_OF_N_GATE
       CASE('EXTERNAL')