Merge remote-tracking branch 'firemodels/master'

Author: cxp484

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -3542,7 +3542,7 @@ \subsection{Louvered Vents}
 \begin{lstlisting}
 &SURF ID='LOUVER', VEL=-2.0, VEL_T=3.0,0.0, TAU_V=5., COLOR='GREEN' /
 \end{lstlisting}
-is a boundary condition for a louvered vent that pushes air into the space with a normal velocity of 2~m/s and a tangential velocity of 3~m/s in the first of the two tangential directions. Note that the negative sign of the normal component of velocity indicates that the fluid is injected into the computational domain. The tangential velocity of 3~m/s indicates that the flow is in the positive $y$ direction. Both the normal and tangential velocity components are ramped up with either \ct{TAU_V} or \ct{RAMP_V}, as shown in Fig.~\ref{fig:tangential_velocity}.
+is a boundary condition for a louvered vent that pushes air into the space with a normal velocity of 2~m/s and a tangential velocity of 3~m/s in the first of the two tangential directions. Note that the negative sign of the normal component of velocity indicates that the fluid is injected into the computational domain. The tangential velocity of 3~m/s indicates that the flow is in the positive $y$ direction. Both the normal and tangential velocity components are ramped up with either \ct{TAU_V} or \ct{RAMP_V}, as shown in Fig.~\ref{fig:tangential_velocity}. If the boundary condition is specified with \ct{MASS_FLUX}, then \ct{VEL_T} should be specified assuming a wall normal velocity of 1~m/s. FDS will then scale the input values for \ct{VEL_T} based on the current wall velocity given by the species mass fractions, mass flux, and wall temperature.
 
 \begin{figure}[ht!]
 \begin{center}
@@ -9140,7 +9140,8 @@ \section{Flux Limiters}
 Central differencing                & \ct{'CENTRAL'}  \\
 Godunov                             & \ct{'GODUNOV'}  \\
 Superbee (VLES, SVLES default)      & \ct{'SUPERBEE'}  \\
-CHARM (DNS, LES default)            & \ct{'CHARM'}  \\
+CHARM (DNS, LES default)            & \ct{'CHARM'}  \\ 
+\hline
 \end{tabular}
 \end{table}
 
@@ -9171,9 +9172,9 @@ \subsection{Density}
 
 The density of the gas has a natural lower bound of zero, but if the density in a cell decreases to nearly zero, the temperature would then increase to an extremely high value due to the equation of state. Thus, by default, the density is kept within the following range:
 \be
-   \min \left( 0.1 \, \rho_\infty , \frac{2 \, p_\infty \, W_{\min}}{\R \, T_{\max}} \right) < \rho <  \frac{2 \, p_\infty \, W_{\max}}{\R \, T_{\min}}
+   0.01 \; \mbox{kg/m$^3$} < \rho < 100 \; \mbox{kg/m$^3$}
 \ee
-where $W_{\min}$ and $W_{\max}$ are the minimum and maximum values of the molecular weight of the tracked gas species in units of g/mol, and $\R$ is the universal gas constant, 8314.5~J/(kmol$\cdot$K). $T_{\min}$ and $T_{\max}$ are described above. To override the limits of density, specify \ct{MINIMUM_DENSITY} and/or \ct{MAXIMUM_DENSITY} on the \ct{CLIP} line in units of kg/m$^3$.
+To override the limits of density, specify \ct{MINIMUM_DENSITY} and/or \ct{MAXIMUM_DENSITY} on the \ct{CLIP} line in units of kg/m$^3$.
 
 Clipping of density and mass fractions violates mass conservation, so it is preferable to avoid clipping if possible.  As discussed in Sec.~\ref{info:Stability_Constraints}, the time step is set to adhere to the CFL constraints of the flow field.  The proper \ct{DT} combined with flux limiters generally avoids the need for clipping.  Beyond this, FDS then employs a mass redistribution scheme, as discussed in FDS Technical Guide~\cite{FDS_Math_Guide}.  If this fails, there is yet one more attempt to avoid clipping---the time step is decreased by 10~\% ($\dt_{\mathrm{new}} = 0.9 \,\dt$) and the scalar transport equations are reiterated.  This process is carried out a maximum of \ct{CLIP_DT_RESTRICTIONS_MAX} times; the default is 5.  In some very extreme circumstances, this loop can drive the time step into numerical instability range ($\dt/\dt_{\mathrm{init}} < \mbox{\ct{LIMITING_DT_RATIO}}$).  You can control the max number of time step restrictions by setting the parameter \ct{CLIP_DT_RESTRICTIONS_MAX} on the \ct{CLIP} line (set to 0 to bypass the algorithm altogether). The number of restrictions (if any) is noted in the \ct{CHID.out} file for a given time step.
 

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

@@ -1522,7 +1522,7 @@ \subsection{1-D Flow in a Tunnel (\texorpdfstring{\ct{tunnel_const_gamma}}{tunne
 \end{eqnarray}
 Note that $\R=8314.5$~J/kmol/K and the sensible enthalpy, $h_{\rm s}=\int_{T_\infty}^T c_p(T') \, \d T'$. The tube is 10~cm long, 1~mm wide and 4~mm tall, and discretized with 1~mm grid cells. A vertical column of cells contains a volumetric heat source that introduces heat at a rate of $\dq'''=2.5133 \times 10^8$~W/m$^3$ within each 1~mm cell. The velocity, density, pressure and temperature are constant starting at the inlet, and at the heat source, the velocity and temperature jump up, and the density and pressure drop down over the span of a single cell and remain steady until the outlet is reached, where the perturbation pressure is set to zero.
 
-For the case where the ratio of specific heats are constant with $\gamma=1.4$ and $W=28.85$ kg/kmol for air, $c_p=\gamma \, (\R/W) /(\gamma-1)=1008.7$~J/kg/K and $\rho c_p T=\gamma \, p_\infty / (\gamma-1)$. In the second equation of (\ref{div1d}), the second two terms on the right hand side are identically zero; thus, the jump in velocity can be computed
+For the case where the ratio of specific heats are constant with $\gamma=1.4$ and $W=28.85$~kg/kmol for air, $c_p=\gamma \, (\R/W) /(\gamma-1)=1008.7$~J/kg/K and $\rho c_p T=\gamma \, p_\infty / (\gamma-1)$. In the second equation of (\ref{div1d}), the second two terms on the right hand side are identically zero; thus, the jump in velocity can be computed
 \be
    \Delta u = \frac{\gamma-1}{\gamma \, p_\infty} \, \dq''' \, \Delta x \approx 0.7087 \; \hbox{m/s}
 \ee
@@ -1554,23 +1554,20 @@ \subsection{1-D Flow in a Tunnel (\texorpdfstring{\ct{tunnel_const_gamma}}{tunne
 \subsection{1-D Flow in a Tunnel (\texorpdfstring{\ct{tunnel_linear_cp}}{tunnel\_linear\_cp})}
 \label{tunnel_linear_cp}
 
-Considering the same layout and the same environmental conditions of the previous case. Air is blown at a constant rate, $u_\infty=1$~m/s, through a 10~cm long, 1~mm wide and 4~mm tall tube. A steady heat source is located in the middle of the tube, composed of a vertical column of cells containing a volumetric heat source that introduces heat at a rate of $\dq'''=2.5133 \times 10^8$~W/m$^3$ within each 1~mm cell . Ignoring wall friction, gravity, radiation transport, thermal conductivity, viscosity, and assuming that the heat source uniformly spans the cross section of the tube, a set of 1-D equations can be derived for the steady-state velocity, $u(x)$, density, $\rho(x)$, perturbation pressure, $\tp(x)$, and temperature, $T(x)$. This equations are the same used in the previous example in (\ref{mom1d}) and (\ref{div1d}).
+Consider the same layout and conditions of the previous case. Air is blown at a constant rate, $u_\infty=1$~m/s, through a 10~cm long, 1~mm wide and 4~mm tall tube. A steady heat source is located in the middle of the tube, composed of a vertical column of cells containing a volumetric heat source that introduces heat at a rate of $\dq'''=2.5133 \times 10^8$~W/m$^3$ within each 1~mm cell. Ignoring wall friction, gravity, radiation transport, thermal conductivity, viscosity, and assuming that the heat source uniformly spans the cross section of the tube, a set of 1-D equations can be derived for the steady-state velocity, $u(x)$, density, $\rho(x)$, perturbation pressure, $\tp(x)$, and temperature, $T(x)$. These equations are the same used in the previous example in (\ref{mom1d}) and (\ref{div1d}).
 
-Remembering that $\R=8314.5$~J/kmol/K and the sensible enthalpy, $h_{\rm s}=\int_{T_\infty}^T c_p(T') \, \d T'$. The velocity, density, pressure and temperature are constant starting at the inlet, until the heat source, where the velocity and temperature jump up, and the density and pressure drop down over the span of a single cell and remain steady until the outlet is reached, where the perturbation pressure is set to zero.
-
-For the case where the specific heat varies as a function of the temperature. In the second equation of (\ref{div1d}), using the first equation of (\ref{mom1d}) and (\ref{div1d}), and some relationships among the variables is possible to arrive to a simpler expression.
+The velocity, density, pressure and temperature are constant starting at the inlet, and at the heat source the velocity and temperature jump up and the density and pressure drop down over the span of a single cell and remain steady until the outlet is reached, where the perturbation pressure is set to zero. In the second equation of (\ref{div1d}), using the first equation of (\ref{mom1d}) and (\ref{div1d}), and some relationships among the variables it is possible to arrive at a simpler expression.
 \begin{eqnarray}
    \rho c_p T \frac{\d u}{\d x}  &=& \dq''' - u \rho \, \frac{\d}{\d x} (h_{\rm s}) - u h_{\rm s} \, \frac{\d}{\d x} (\rho) - c_p T u \frac{\d \rho}{\d x}  + h_{\rm s} u \frac{\d \rho}{\d x} \\
    c_p \rho_\infty u_\infty \frac{\d T}{\d x} &=& \dq'''  - u \rho \, \frac{\d}{\d x} (h_{\rm s}) + c_p u \rho \frac{\d T}{\d x}  \\
   \dq''' &=&   u \rho \, \frac{\d}{\d x} (h_{\rm s}) \\
   \frac{ \dq'''}{\rho_\infty u_\infty} &=&  \frac{\d}{\d x} (h_{\rm s})  \label{1stLaw}
 \end{eqnarray}
-
-Through the steps, the continuity, in (\ref{mom1d}), and state equation, in (\ref{div1d}), are used to make the temperature the only variable and simplifying equivalent terms until the equation (\ref{1stLaw}) is obtained. Equation (\ref{1stLaw}) corresponds to a 1-D form of the first law of thermodynamics, proving the divergence equation as a consistent energy balance. Then, for the case where the $c_p$ is a linear function of the temperature, $c_p=aT +b$, we have
+The continuity, (\ref{mom1d}), and state equation, (\ref{div1d}), render the temperature the only variable and after simplifying, Eq.~(\ref{1stLaw}) is obtained. Equation (\ref{1stLaw}) corresponds to a 1-D form of the first law of thermodynamics, proving the divergence equation as a consistent energy balance. Then, for the case where the $c_p$ is a linear function of the temperature, $c_p=aT +b$, we have
 \begin{eqnarray}
   \frac{ \dq'''}{\rho_\infty u_\infty} \Delta x &=&  T^2 \frac{a}{2} +bT -  (T_\infty^2 \frac{a}{2} +bT_\infty) \label{LinearCP}
 \end{eqnarray}
- Substituting in Eq.~(\ref{LinearCP}) for $a=0.1584$ J/kg/K$^2$ and $b=953.5650$ J/kg/K (these coefficients give $c_p(T_\infty)=1000$ J/kg/K) the temperature downstream of the heat source is computed to be 499.2~K. Correspondingly, the density and velocity are, $\rho_1$= 0.7046~\si{kg/m^3} and $u_1= 1.7030$~m/s. To determine the drop in pressure, Eq.~(\ref{mom1d}) is written in discretized form as in the previous verification. Then assuming that the subscript 0 refers to the cell containing the heat source, 1 to the first cell downstream, and -1 to the first cell upstream. The velocity in cell 0 is the approximated as the average of the velocity upstream and downstream, $u_0=1.3515$~m/s.
+Substituting in Eq.~(\ref{LinearCP}) for $a=0.1584$~J/kg/K$^2$ and $b=953.5650$~J/kg/K (these coefficients give $c_p(T_\infty)=1000$~J/kg/K) the temperature downstream of the heat source is computed to be 499.2~K. Correspondingly, the density and velocity are, $\rho_1=0.7046$~\si{kg/m^3} and $u_1= 1.7030$~m/s. To determine the drop in pressure, Eq.~(\ref{mom1d}) is written in discretized form as in the previous verification. Then assuming that the subscript 0 refers to the cell containing the heat source, 1 to the first cell downstream, and -1 to the first cell upstream. The velocity in cell 0 is the approximated as the average of the velocity upstream and downstream, $u_0=1.3515$~m/s.
 \begin{eqnarray}
    \tp_1 - \tp_0    &=& - \frac{u_1^2 - u_0^2}{2}   \rho_1    \approx -0.3783 \; \hbox{Pa} \\
    \tp_0 - \tp_{-1} &=& - \frac{u_0^2 - u_{-1}^2}{2}  \rho_0  \approx -0.3936 \; \hbox{Pa}
@@ -1669,7 +1666,7 @@ \subsection{Dynamic Smagorinsky (\texorpdfstring{\ct{dsmag}}{dsmag})}
 
 In the previous section, all calculations were performed with a constant and uniform Smagorinsky coefficient, $C_{\rm s} = 0.2$.  For the canonical case of homogeneous decaying isotropic turbulence -- \emph{at sufficiently high Reynolds number} -- this model is sufficient.  However, we noticed that even for the isotropic turbulence problem when the grid Reynolds number is low (i.e., the flow is well-resolved) the constant coefficient model tends to over predict the dissipation of kinetic energy (see Fig.~\ref{fig_cbc_energy}).  This is because the eddy viscosity does not converge to zero at the proper rate; so long as strain is present in the flow (the magnitude of the strain rate tensor is nonzero), the eddy viscosity will be nonzero.  This violates a guiding principle in LES development: that the method should converge to a DNS if the flow field is sufficiently resolved.
 
-The dynamic procedure for calculating the model coefficient (set \ct{TURBULENCE_MODEL=`DYNAMIC SMAGORINSKY'} on \ct{MISC}) alleviates this problem.  The basis of the model is that the coefficient should be the same for two different filter scales within the inertial subrange.  Details of the procedure are explained in the following references \cite{Germano:1991,PinoMartin:2000,Moin:1991,Lund:1997,FDS_Math_Guide}.  Here we present results for the implementation of the dynamic model in FDS.  In Fig.~\ref{fig_c_smag} we show contours of the Smagorinsky coefficient $C_{\rm s}(\mathbf{x},t)$ at a time midway through a $64^3$ simulation of the CBC experiment.  Notice that the coefficient ranges from 0.00 to roughly 0.30 within the domain with the average value falling around 0.17.
+The dynamic procedure for calculating the model coefficient (set \ct{TURBULENCE_MODEL='DYNAMIC SMAGORINSKY'} on \ct{MISC}) alleviates this problem.  The basis of the model is that the coefficient should be the same for two different filter scales within the inertial subrange.  Details of the procedure are explained in the following references \cite{Germano:1991,PinoMartin:2000,Moin:1991,Lund:1997,FDS_Math_Guide}.  Here we present results for the implementation of the dynamic model in FDS.  In Fig.~\ref{fig_c_smag} we show contours of the Smagorinsky coefficient $C_{\rm s}(\mathbf{x},t)$ at a time midway through a $64^3$ simulation of the CBC experiment.  Notice that the coefficient ranges from 0.00 to roughly 0.30 within the domain with the average value falling around 0.17.
 
 \begin{figure}[t]
 \centering
@@ -2392,7 +2389,7 @@ \subsubsection{Case 4: Gases Released by Solid Phase Reactions at the Boundary o
 \subsection{Mass Flux through Domain Boundaries (\texorpdfstring{\ct{mass_flux_wall}}{mass\_flux\_wall})}
 \label{mass_flux_wall}
 
-This set of tests considers a specified mass injection of water vapor from a \ct{VENT} on the floor of a compartment.  The entire top of the compartment is \ct{OPEN} to the ambient.  The lighter water vapors are buoyant, so they rise creating a turbulent flow into and out of the domain.  Therefore, the input is steady, while the outflow of water vapor is quite chaotic.  To judge whether a proper mass balance is achieved by FDS we must monitor not only the inflow and outflow, but also the total mass of water vapor in the domain at a given time.  Therefore, we add \ct{DEVC} of quantity \ct{MASS_FLUX_WALL} and \ct{STATISTICS=`SURFACE INTEGRAL`} for inflow and outflow, and we also add \ct{MASS_FILE=.TRUE.} to the \ct{DUMP} line to monitor the total water vapor at any point in time in the domain.  With $\rho$ and $Y_{\mbox{\tiny H2O}}$ denoting the local mass density and mass fraction of water vapor in a cell, and with $\dot{m}_{\mbox{\tiny H2O}}^{\prime\prime}$ denoting the flux of water vapor at a cell face on the domain boundary, the mass balance may be written as,
+This set of tests considers a specified mass injection of water vapor from a \ct{VENT} on the floor of a compartment.  The entire top of the compartment is \ct{OPEN} to the ambient.  The lighter water vapors are buoyant, so they rise creating a turbulent flow into and out of the domain.  Therefore, the input is steady, while the outflow of water vapor is quite chaotic.  To judge whether a proper mass balance is achieved by FDS we must monitor not only the inflow and outflow, but also the total mass of water vapor in the domain at a given time.  Therefore, we add \ct{DEVC} of quantity \ct{MASS_FLUX_WALL} and \ct{STATISTICS='SURFACE INTEGRAL'} for inflow and outflow, and we also add \ct{MASS_FILE=.TRUE.} to the \ct{DUMP} line to monitor the total water vapor at any point in time in the domain.  With $\rho$ and $Y_{\mbox{\tiny H2O}}$ denoting the local mass density and mass fraction of water vapor in a cell, and with $\dot{m}_{\mbox{\tiny H2O}}^{\prime\prime}$ denoting the flux of water vapor at a cell face on the domain boundary, the mass balance may be written as,
 \begin{equation}
 \label{eq:H2O_mass_balance}
 \frac{\mbox{d}}{\mbox{d}t} \left( \int_{V} \rho Y_{\mbox{\tiny H2O}} \,\mbox{d}V \right) + \int_{S_{out}} \dot{m}_{\mbox{\tiny H2O}}^{\prime\prime} \,\mbox{d}S - \int_{S_{in}} \dot{m}_{\mbox{\tiny H2O}}^{\prime\prime} \,\mbox{d}S = 0
@@ -2412,7 +2409,7 @@ \subsection{Mass Flux through Domain Boundaries (\texorpdfstring{\ct{mass_flux_w
 \subsection{Mass Balance on a Gas Control Volume (\texorpdfstring{\ct{mass_balance_gas_volume}}{mass\_balance\_gas\_volume})}
 \label{mass_balance_gas_volume}
 
-Using the same computational setup as the \ct{mass_flux_wall} series above, the \ct{mass_balance_gas_volume} test case defines a control volume (CV) internal to the computational domain, 6 m on a side (from 2 m to 8 m in each direction).  The mass of water in the CV is output using a \ct{SPATIAL_STATISTIC=`VOLUME INTEGRAL`} with \ct{QUANTITY=`DENSITY`} of the lumped species for water and the flows into and out of the CV are obtained using \ct{QUANTITY=`TOTAL MASS FLUX X`}, etc., with \ct{SPATIAL_STATISTIC=`AREA INTEGRAL`}. Figure \ref{fig:mass_balance_gas_volume} shows the resulting balance of accumulation (dm/dt) and the inflow and outflow of the control volume.
+Using the same computational setup as the \ct{mass_flux_wall} series above, the \ct{mass_balance_gas_volume} test case defines a control volume (CV) internal to the computational domain, 6 m on a side (from 2 m to 8 m in each direction).  The mass of water in the CV is output using a \ct{SPATIAL_STATISTIC='VOLUME INTEGRAL'} with \ct{QUANTITY='DENSITY'} of the lumped species for water and the flows into and out of the CV are obtained using \ct{QUANTITY='TOTAL MASS FLUX X'}, etc., with \ct{SPATIAL_STATISTIC='AREA INTEGRAL'}. Figure \ref{fig:mass_balance_gas_volume} shows the resulting balance of accumulation (dm/dt) and the inflow and outflow of the control volume.
 
 \begin{figure}[ht]
 \centering

Source/data.f90

@@ -2452,7 +2452,7 @@ SUBROUTINE THERMO_TABLE_LIQUID(I_TMP,CP,H_V,T_REF,T_MELT,T_BOIL,SPEC_INDEX,FUEL,
 LOGICAL,INTENT(OUT) :: FUEL
 TYPE (THERMO_DATA_TYPE), POINTER :: TD
 
-!Default to water properteis
+! Default to water properties
 
 IF (SPEC_INDEX < 0) THEN
    FUEL = .FALSE.

Source/vege.f90

@@ -430,9 +430,9 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
                U_Z(2) = 0.5_EB*(U(IIG-1,JJG,KWIND)+U(IIG,JJG,KWIND))
                V_Z(1) = 0.5_EB*(V(IIG,JJG-1,KWIND-1)+V(IIG,JJG,KWIND-1))
                V_Z(2) = 0.5_EB*(V(IIG,JJG-1,KWIND)+V(IIG,JJG,KWIND))
-               ! If wind comes from first grid cell assume zero wind at ground level
+               ! If wind comes from first grid cell assume plug flow near ground level
                IF (KWIND==1) THEN
-                  U_Z(1) = 0._EB; V_Z(1) = 0._EB; ZWIND(1) = 0._EB
+                  U_Z(1) = U_Z(2); V_Z(1) = V_Z(2); ZWIND(1) = Z_LS(IIG,JJG)
                ENDIF
                U_LS(IIG,JJG) = U_Z(1) + (REF_WIND_HEIGHT-ZWIND(1))/(ZWIND(2)-ZWIND(1))*(U_Z(2)-U_Z(1))
                V_LS(IIG,JJG) = V_Z(1) + (REF_WIND_HEIGHT-ZWIND(1))/(ZWIND(2)-ZWIND(1))*(V_Z(2)-V_Z(1))