Skip to content

Commit bf7682b

Browse files
committed
Merge remote-tracking branch 'firemodels/master' into FireX
2 parents 5ab8b47 + 65fa1c9 commit bf7682b

File tree

16 files changed

+724
-41
lines changed

16 files changed

+724
-41
lines changed

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

Lines changed: 82 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -5025,24 +5025,99 @@ \section{SFPE Verification Cases}
50255025
\subsection{Case 6: 2-D Heat Transfer with Cooling by Convection}
50265026
\label{SFPE_Case_6}
50275027

5028-
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.
5028+
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}).
50295029

5030-
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}.
5030+
\begin{figure}[ht]
5031+
\centering
5032+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_6}
5033+
\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.}
5034+
\label{fig:SFPE_Case_6}
5035+
\end{figure}
5036+
5037+
\FloatBarrier
50315038

50325039
\subsection{Case 7: 2-D Heat Transfer by Convection and Radiation}
50335040
\label{SFPE_Case_7}
50345041

50355042
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
50365043
\be
5037-
T(t) = T_\infty + 345 \, \ln \left( 8t/60 + 1 \right)
5044+
T(t) = T_\infty + 345 \, \ln \left( 8t/60 + 1 \right) \label{ISO_834}
50385045
\ee
5039-
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}.
5046+
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}).
50405047

50415048
\begin{figure}[ht]
5042-
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_6}
5049+
\centering
50435050
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_7}
5044-
\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.}
5045-
\label{fig:SFPE_Case_6_7}
5051+
\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.}
5052+
\label{fig:SFPE_Case_7}
5053+
\end{figure}
5054+
5055+
\FloatBarrier
5056+
5057+
5058+
\subsection{Case 8: 2-D Heat Transfer with Temperature-Dependent Conductivity}
5059+
\label{SFPE_Case_8}
5060+
5061+
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}).
5062+
5063+
\begin{figure}[ht]
5064+
\centering
5065+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_8}
5066+
\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.}
5067+
\label{fig:SFPE_Case_8}
5068+
\end{figure}
5069+
5070+
\FloatBarrier
5071+
5072+
5073+
\subsection{Case 9: 2-D Heat Transfer in a Composite Section with Temperature-Dependent Conductivity}
5074+
\label{SFPE_Case_9}
5075+
5076+
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}).
5077+
5078+
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.
5079+
5080+
\begin{figure}[ht]
5081+
\centering
5082+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_9}
5083+
\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.}
5084+
\label{fig:SFPE_Case_9}
5085+
\end{figure}
5086+
5087+
\FloatBarrier
5088+
5089+
\subsection{Case 13: 2-D Heat Transfer with Moisture Evaporation}
5090+
\label{SFPE_Case_13}
5091+
5092+
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}).
5093+
5094+
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
5095+
\be
5096+
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
5097+
\ee
5098+
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.
5099+
5100+
\begin{figure}[ht]
5101+
\centering
5102+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_13b}
5103+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_13}
5104+
\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.}
5105+
\label{fig:SFPE_Case_13}
5106+
\end{figure}
5107+
5108+
5109+
\FloatBarrier
5110+
5111+
\subsection{Case 16: 3-D Heat Transfer with Non-Uniform Heat Flux}
5112+
\label{SFPE_Case_16}
5113+
5114+
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}).
5115+
5116+
\begin{figure}[ht]
5117+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_16a}
5118+
\includegraphics[height=2.2in]{SCRIPT_FIGURES/SFPE_Case_16b}
5119+
\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.}
5120+
\label{fig:SFPE_Case_16}
50465121
\end{figure}
50475122

50485123

Source/data.f90

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1280,21 +1280,21 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
12801280
OUTPUT_QUANTITY(534)%SHORT_NAME = 'rho*(Z*u + D*dZ/dx)'
12811281
OUTPUT_QUANTITY(534)%CELL_POSITION = CELL_FACE
12821282
OUTPUT_QUANTITY(534)%IOR = 1
1283-
OUTPUT_QUANTITY(534)%SPEC_ID_REQUIRED=.TRUE.
1283+
!OUTPUT_QUANTITY(534)%SPEC_ID_REQUIRED=.TRUE.
12841284

12851285
OUTPUT_QUANTITY(535)%NAME = 'TOTAL MASS FLUX Y'
12861286
OUTPUT_QUANTITY(535)%UNITS = 'kg/s/m2'
12871287
OUTPUT_QUANTITY(535)%SHORT_NAME = 'rho*(Z*v + D*dZ/dy)'
12881288
OUTPUT_QUANTITY(535)%CELL_POSITION = CELL_FACE
12891289
OUTPUT_QUANTITY(535)%IOR = 2
1290-
OUTPUT_QUANTITY(535)%SPEC_ID_REQUIRED=.TRUE.
1290+
!OUTPUT_QUANTITY(535)%SPEC_ID_REQUIRED=.TRUE.
12911291

12921292
OUTPUT_QUANTITY(536)%NAME = 'TOTAL MASS FLUX Z'
12931293
OUTPUT_QUANTITY(536)%UNITS = 'kg/s/m2'
12941294
OUTPUT_QUANTITY(536)%SHORT_NAME = 'rho*(Z*w + D*dZ/dz)'
12951295
OUTPUT_QUANTITY(536)%CELL_POSITION = CELL_FACE
12961296
OUTPUT_QUANTITY(536)%IOR = 3
1297-
OUTPUT_QUANTITY(536)%SPEC_ID_REQUIRED=.TRUE.
1297+
!OUTPUT_QUANTITY(536)%SPEC_ID_REQUIRED=.TRUE.
12981298

12991299
OUTPUT_QUANTITY(550)%NAME = 'CUTCELL VELOCITY DIVERGENCE'
13001300
OUTPUT_QUANTITY(550)%UNITS = '1/s'

Source/dump.f90

Lines changed: 27 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -9657,29 +9657,40 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
96579657
ENDIF
96589658
CASE(534) ! TOTAL MASS FLUX X
96599659
GAS_PHASE_OUTPUT_RES = 0._EB
9660-
IF (Z_INDEX>0) THEN
9661-
GAS_PHASE_OUTPUT_RES = ADV_FX(II,JJ,KK,Z_INDEX) + DIF_FX(II,JJ,KK,Z_INDEX)
9662-
ELSEIF (Y_INDEX>0) THEN
9663-
GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
9664-
(ADV_FX(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FX(II,JJ,KK,1:N_TRACKED_SPECIES)) )
9660+
IF (Z_INDEX<=0 .AND. Y_INDEX<=0) THEN
9661+
GAS_PHASE_OUTPUT_RES = SUM(ADV_FX(II,JJ,KK,:) + DIF_FX(II,JJ,KK,:))
9662+
ELSE
9663+
IF (Z_INDEX>0) THEN
9664+
GAS_PHASE_OUTPUT_RES = ADV_FX(II,JJ,KK,Z_INDEX) + DIF_FX(II,JJ,KK,Z_INDEX)
9665+
ELSEIF (Y_INDEX>0) THEN
9666+
GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
9667+
(ADV_FX(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FX(II,JJ,KK,1:N_TRACKED_SPECIES)) )
9668+
ENDIF
96659669
ENDIF
96669670
CASE(535) ! TOTAL MASS FLUX Y
96679671
GAS_PHASE_OUTPUT_RES = 0._EB
9668-
IF (Z_INDEX>0) THEN
9669-
GAS_PHASE_OUTPUT_RES = ADV_FY(II,JJ,KK,Z_INDEX) + DIF_FY(II,JJ,KK,Z_INDEX)
9670-
ELSEIF (Y_INDEX>0) THEN
9671-
GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
9672-
(ADV_FY(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FY(II,JJ,KK,1:N_TRACKED_SPECIES)) )
9672+
IF (Z_INDEX<=0 .AND. Y_INDEX<=0) THEN
9673+
GAS_PHASE_OUTPUT_RES = SUM(ADV_FY(II,JJ,KK,:) + DIF_FY(II,JJ,KK,:))
9674+
ELSE
9675+
IF (Z_INDEX>0) THEN
9676+
GAS_PHASE_OUTPUT_RES = ADV_FY(II,JJ,KK,Z_INDEX) + DIF_FY(II,JJ,KK,Z_INDEX)
9677+
ELSEIF (Y_INDEX>0) THEN
9678+
GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
9679+
(ADV_FY(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FY(II,JJ,KK,1:N_TRACKED_SPECIES)) )
9680+
ENDIF
96739681
ENDIF
96749682
CASE(536) ! TOTAL MASS FLUX Z
96759683
GAS_PHASE_OUTPUT_RES = 0._EB
9676-
IF (Z_INDEX>0) THEN
9677-
GAS_PHASE_OUTPUT_RES = ADV_FZ(II,JJ,KK,Z_INDEX) + DIF_FZ(II,JJ,KK,Z_INDEX)
9678-
ELSEIF (Y_INDEX>0) THEN
9679-
GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
9680-
(ADV_FZ(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FZ(II,JJ,KK,1:N_TRACKED_SPECIES)) )
9684+
IF (Z_INDEX<=0 .AND. Y_INDEX<=0) THEN
9685+
GAS_PHASE_OUTPUT_RES = SUM(ADV_FZ(II,JJ,KK,:) + DIF_FZ(II,JJ,KK,:))
9686+
ELSE
9687+
IF (Z_INDEX>0) THEN
9688+
GAS_PHASE_OUTPUT_RES = ADV_FZ(II,JJ,KK,Z_INDEX) + DIF_FZ(II,JJ,KK,Z_INDEX)
9689+
ELSEIF (Y_INDEX>0) THEN
9690+
GAS_PHASE_OUTPUT_RES = DOT_PRODUCT( Z2Y(Y_INDEX,1:N_TRACKED_SPECIES),&
9691+
(ADV_FZ(II,JJ,KK,1:N_TRACKED_SPECIES) + DIF_FZ(II,JJ,KK,1:N_TRACKED_SPECIES)) )
9692+
ENDIF
96819693
ENDIF
9682-
96839694
CASE(550) ! CUTCELL VELOCITY DIVERGENCE
96849695
GAS_PHASE_OUTPUT_RES = CARTVELDIV(II,JJ,KK)
96859696

Source/pres.f90

Lines changed: 15 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -2896,17 +2896,23 @@ SUBROUTINE ULMAT_H_MATRIX_SOLVER_SETUP(NM,IPZ)
28962896
CALL HYPRE_IJMATRIXSETOBJECTTYPE(ZM%HYPRE_ZM%A_H, HYPRE_PARCSR, HYPRE_IERR)
28972897
CALL HYPRE_IJMATRIXINITIALIZE_V2(ZM%HYPRE_ZM%A_H, HYPRE_MEMORY_HOST, HYPRE_IERR)
28982898
IF (ZM%MTYPE==SYMM_INDEFINITE) THEN
2899-
! Rows 1 to ZM%NUNKH-1, last column, set all to zero:
2900-
DO IROW=1,ZM%NUNKH-1
2899+
IF(ZM%NUNKH==1) THEN ! Single unknown, zero coefficient matrix (1,1). Set coefficient to 1.
2900+
NNZ_H_MAT(1) = 1
2901+
DEALLOCATE(JD_H_MAT); ALLOCATE(JD_H_MAT(1,1)); JD_H_MAT(1,1) = 1
2902+
DEALLOCATE( D_H_MAT); ALLOCATE( D_H_MAT(1,1)); D_H_MAT(1,1) = 1._EB
2903+
ELSE ! More than one unknown
2904+
! Rows 1 to ZM%NUNKH-1, last column, set all to zero:
2905+
DO IROW=1,ZM%NUNKH-1
2906+
DO JCOL=1,NNZ_H_MAT(IROW)
2907+
IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) CYCLE ! Make zero matrix entries in last column.
2908+
D_H_MAT(JCOL,IROW) = 0._EB
2909+
ENDDO
2910+
ENDDO
2911+
! Last row, all zeros except the diagonal that keeps diagonal number: Note after previous loop IROW==ZM%NUNKH
29012912
DO JCOL=1,NNZ_H_MAT(IROW)
2902-
IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) CYCLE ! Make zero matrix entries in last column.
2903-
D_H_MAT(JCOL,IROW) = 0._EB
2913+
IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) D_H_MAT(JCOL,IROW) = 0._EB
29042914
ENDDO
2905-
ENDDO
2906-
! Last row, all zeros except the diagonal that keeps diagonal number: Note after previous loop IROW==ZM%NUNKH
2907-
DO JCOL=1,NNZ_H_MAT(IROW)
2908-
IF ( JD_H_MAT(JCOL,IROW) /= ZM%NUNKH ) D_H_MAT(JCOL,IROW) = 0._EB
2909-
ENDDO
2915+
ENDIF
29102916
ENDIF
29112917
JD_H_MAT = JD_H_MAT - 1
29122918
DO IROW=1,ZM%NUNKH

Source/read.f90

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14515,7 +14515,6 @@ SUBROUTINE READ_CTRL
1451514515
CF%CONTROL_INDEX = OR_GATE
1451614516
CASE('ONLY')
1451714517
CF%CONTROL_INDEX = XOR_GATE
14518-
CF%N = 1
1451914518
CASE('AT_LEAST')
1452014519
CF%CONTROL_INDEX = X_OF_N_GATE
1452114520
CASE('EXTERNAL')

0 commit comments

Comments
 (0)