Merge pull request #16423 from cxp484/FireX

FireX: Merge with firemodels/master

Author: cxp484

Manuals/FDS_Technical_Reference_Guide/FIGURES/Rad_rand_rot_cube_hit.pdf

@@ -500,10 +500,13 @@ \subsection{Spatial Discretization}
 where the $D$-terms on the LHS are integrals (\ref{cellfaceintegral}) evaluated on upwind and downwind sides of the cell. In the rectilinear mesh,
 $D_{xu}^l = -D_{xd}^l$ and the equation can be simplified further. In addition, the integrals over the solid angle can be calculated analytically
 \begin{align}
+\label{eq:directional_weights1}
    D_x^l & = \int_{\Omega^l}(\bs^l\cdot {\bf i}) \, \d\Omega = \int_{\dph}\int_{\dth} (\bs^l\cdot{\bf i}) \sin\theta \; \d \theta \d\phi = \int_{\dph}\int_{\dth} \cos\phi\sin\theta\sin\theta\; \d\theta \d\phi \notag \\
          & = \frac{1}{2}\left( \sin\phi^+ - \sin\phi^- \right) \left[\Delta\theta - \left(\cos\theta^+\sin\theta^+ - \cos\theta^-\sin\theta^-\right)\right]  \\
+\label{eq:directional_weights2}
    D_y^l & = \int_{\Omega^l}(\bs^l\cdot {\bf j}) \, \d\Omega  =  \int_{\dph}\int_{\dth} \sin\phi\sin\theta\sin\theta\; \d\theta \d\phi \notag \\
          & = \frac{1}{2}\left( \cos\phi^- - \cos\phi^+ \right) \left[\Delta\theta - \left(\cos\theta^+\sin\theta^+ - \cos\theta^-\sin\theta^-\right)\right]  \\
+\label{eq:directional_weights3}
    D_z^l & = \int_{\Omega^l}(\bs^l\cdot {\bf k}) \d \Omega = \int_{\dph}\int_{\dth} \cos\theta\sin\theta\;  \d\theta \d\phi \notag \\
          & = \frac{1}{2}\Delta\phi \left[\left(\sin\theta^+\right)^2 - \left(\sin\theta^-\right)^2\right]
 \\
@@ -535,7 +538,225 @@ \subsection{Spatial Discretization}
 all the upwind intensities are determined from solid or gas phase boundaries. In theory, iterations are needed if the reflections
 or scattering are important, or if the scenario is optically very thick. Currently, no iterations are made.
 
-To mitigate ray effects in the wall radiative heat flux, an option is implemented to randomly rotate the local meridian ($x$), azimuthal ($y$), and axial ($z$) coordinate system at each selected time step with respect to which the angular coordinates, $\theta$ and $\phi$ are evaluated. When random rotation is enabled, the discrete angular integration coefficients (e.g., $D_{xu}^l$) are no longer exactly symmetric with respect to the mesh faces. Consequently, the sum of $D_{xu}^l$ (and the corresponding coefficients for the other coordinate directions) over a face may deviate slightly from the theoretical value of $\pi$.
+\textbf{Random rotation of angular directions:}
+To mitigate ray effects in the predicted wall radiative heat flux, the local meridional ($x'$), azimuthal ($y'$), and axial ($z'$) coordinate system may be randomly rotated at selected time steps. The angular coordinates, $\theta$ and $\phi$, are then evaluated with respect to the rotated coordinate system. Figure~\ref{fig:FVMAngularRotation} illustrates the angular discretization for several orientations of the local coordinate system. In this example, the azimuthal ($y'$) and axial ($z'$) axes are rotated in the $YZ$ plane, while the meridional ($x'$) axis remains unchanged. When random rotation is enabled, a new random orientation of the local coordinate system is assigned at each selected time step. Consequently, rays originating from a source intersect distant target surfaces along different discrete directions over time. This variation improves the effective angular coverage and substantially reduces ray effects in the time-averaged wall radiative heat flux. Figure~\ref{fig:rad_rand_rot_cube_hit} demonstrates this effect for an angular discretization containing 24 discrete directions. After 100 random rotations, the ray-intersection locations provide nearly uniform coverage of the enclosure surfaces. Increasing the number of rotations to 500 produces still more uniform surface coverage.
+
+With this rotation, the angular coordinates $\theta$ and $\phi$ are defined with respect to the rotated coordinate system $(x',y',z')$. Let the unit vectors of the rotated meridional, azimuthal, and axial axes be denoted by $\mathbf{v}_{\mathrm{mer}}$, $\mathbf{v}_{\mathrm{azi}}$, and $\mathbf{v}_{\mathrm{axi}}$, respectively. Then, for the $l$th discrete direction with local direction cosines $(\mu_{\mathrm{mer},l},\mu_{\mathrm{azi},l},\mu_{\mathrm{axi},l})$, the corresponding direction vector in the global Cartesian coordinate system is given by
+
+\begin{equation}
+\begin{aligned}
+\begin{bmatrix}
+D_x^l\\
+D_y^l\\
+D_z^l
+\end{bmatrix}
+&=
+\mu_{\mathrm{mer},l}\,\mathbf{v}_{\mathrm{mer}}
++\mu_{\mathrm{azi},l}\,\mathbf{v}_{\mathrm{azi}}
++\mu_{\mathrm{axi},l}\,\mathbf{v}_{\mathrm{axi}} \\
+&=
+\begin{bmatrix}
+v_{\mathrm{mer},x} & v_{\mathrm{azi},x} & v_{\mathrm{axi},x}\\
+v_{\mathrm{mer},y} & v_{\mathrm{azi},y} & v_{\mathrm{axi},y}\\
+v_{\mathrm{mer},z} & v_{\mathrm{azi},z} & v_{\mathrm{axi},z}
+\end{bmatrix}
+\begin{bmatrix}
+\mu_{\mathrm{mer},l}\\
+\mu_{\mathrm{azi},l}\\
+\mu_{\mathrm{axi},l}
+\end{bmatrix},
+\end{aligned}
+\label{eq:rotated_direction_cosines}
+\end{equation}
+
+where the local direction cosines $(\mu_{\mathrm{mer},l},\mu_{\mathrm{azi},l},\mu_{\mathrm{axi},l})$ are obtained from Eqs.~\ref{eq:directional_weights1}--\ref{eq:directional_weights3}.
+
+When rotation is enabled, the discrete angular integration coefficients (e.g., $D_{xu}^l$) are no longer exactly symmetric with respect to the mesh faces. Consequently, the summation of $D_{xu}^l$ (and the corresponding coefficients for the other coordinate directions) over a face deviates slightly from the theoretical value of $\pi$. For the default discretization of 100 solid angles, these deviations are typically very small and decrease further as the angular resolution is increased.
+
+\newcommand{\DrawRADFVMSphere}[1]{%
+   \pgfmathsetmacro{\ang}{#1}
+
+   % Sphere
+   \shade[ball color=blue!10] (0,0) circle (1);
+   \draw[thick] (0,0) circle (1);
+
+   % Latitude lines
+   \foreach \i in {0,...,\NRT}
+   {
+     \pgfmathsetmacro{\theta}{180*\i/\NRT}
+
+     \draw[blue!60, thick]
+     plot[
+         variable=\phi,
+         domain=0:360,
+         samples=150,
+         smooth
+     ]
+     (
+     {
+         (sin(\theta)*sin(\phi))*cos(\ang)
+         +(cos(\theta))*sin(\ang)
+     },
+     {
+        -(sin(\theta)*sin(\phi))*sin(\ang)
+         +(cos(\theta))*cos(\ang)
+     }
+     );
+   }
+
+   % Longitude lines
+   \foreach \i in {0,...,7}
+   {
+     \pgfmathsetmacro{\thetaLow}{180*\i/\NRT}
+     \pgfmathsetmacro{\thetaHigh}{180*(\i+1)/\NRT}
+
+     \pgfmathsetmacro{\nphiRaw}{
+         0.5*\NRA*(cos(\thetaLow)-cos(\thetaHigh))
+     }
+
+     \pgfmathtruncatemacro{\nphi}{
+         max(4,4*round(\nphiRaw/4))
+     }
+
+     \foreach \j in {0,...,\numexpr\nphi-1\relax}
+     {
+         \pgfmathsetmacro{\phi}{360*\j/\nphi}
+
+         \draw[red!60, thick]
+         plot[
+             variable=\theta,
+             domain=\thetaLow:\thetaHigh,
+             samples=30,
+             smooth
+         ]
+         (
+         {
+             (sin(\theta)*sin(\phi))*cos(\ang)
+             +(cos(\theta))*sin(\ang)
+         },
+         {
+            -(sin(\theta)*sin(\phi))*sin(\ang)
+             +(cos(\theta))*cos(\ang)
+         }
+         );
+     }
+   }
+
+   %-----------------------------------
+   % Axes
+   %-----------------------------------
+
+   \ifnum\ang=0
+
+       % No rotation: y=y', z=z'
+       \draw[->, very thick, black]
+           (0,0) -- (1.45,0)
+           node[right] {$y,y'$};
+
+       \draw[->, very thick, black]
+           (0,0) -- (0,1.45)
+           node[above] {$z,z'$};
+
+   \else
+
+       %-----------------------------------
+       % Default fixed axes
+       %-----------------------------------
+
+       \draw[->, thick, black]
+           (0,0) -- (1.45,0)
+           node[right] {$y$};
+
+       \draw[->, thick, black]
+           (0,0) -- (0,1.45)
+           node[above] {$z$};
+
+
+       %-----------------------------------
+       % Rotated axes
+       %-----------------------------------
+
+       \draw[->, very thick, red]
+           (0,0) --
+           ({1.45*cos(\ang)},
+            {-1.45*sin(\ang)})
+           node[pos=1.05, right] {$y'$};
+
+       \draw[->, very thick, red]
+           (0,0) --
+           ({1.45*sin(\ang)},
+            {1.45*cos(\ang)})
+           node[pos=1.05, above] {$z'$};
+
+   \fi
+
+
+   %-----------------------------------
+   % X-axis (out of page)
+   %-----------------------------------
+
+   \draw[thick] (0,0) circle (0.05);
+   \fill (0,0) circle (0.015);
+   \node[below left] at (0,0) {$x,x'$};
+
+}
+
+\begin{figure}[htbp]
+\centering
+
+\begin{tikzpicture}[scale=2.3]
+
+\def\NRA{100}
+\pgfmathtruncatemacro{\NRT}{
+    2*round(0.5*1.17*\NRA^(1/2.26))
+}
+
+%-------------------------
+% Top left
+%-------------------------
+\begin{scope}[shift={(0,4)}]
+    \DrawRADFVMSphere{0}
+    \node at (0,-1.65) {(a) $0^\circ$};
+\end{scope}
+
+%-------------------------
+% Top right
+%-------------------------
+\begin{scope}[shift={(3.5,4)}]
+    \DrawRADFVMSphere{45}
+    \node at (0,-1.65) {(b) $45^\circ$};
+\end{scope}
+
+%-------------------------
+% Bottom left
+%-------------------------
+\begin{scope}[shift={(0,0)}]
+    \DrawRADFVMSphere{135}
+    \node at (0,-1.65) {(c) $135^\circ$};
+\end{scope}
+
+%-------------------------
+% Bottom right
+%-------------------------
+\begin{scope}[shift={(3.5,0)}]
+    \DrawRADFVMSphere{225}
+    \node at (0,-1.65) {(d) $225^\circ$};
+\end{scope}
+
+\end{tikzpicture}
+
+\caption[Finite volume method angular discretization]{Finite volume method angular discretization with the default 100 discrete directions (which turned out to be 104). The angular discretization consists of 8 divisions in the polar angle $\theta$, with the corresponding number of divisions in the azimuthal angle $\phi$ given by [4, 12, 16, 20, 20, 16, 12, 4]. Top left: default angular discretization without rotation. Top right: $45^\circ$ clockwise rotation of the azimuthal ($y'$) and axial ($z'$) coordinate axes in the $YZ$ plane, relative to which the angular discretization is performed. Bottom left: $125^\circ$ clockwise rotation. Bottom right: $225^\circ$ clockwise rotation.}
+\label{fig:FVMAngularRotation}
+\end{figure}
+
+\begin{figure}[ht]
+\begin{center}
+\includegraphics[height=5in]{FIGURES/Rad_rand_rot_cube_hit}
+\caption[Random rotation of radiation angles]{Distribution of radiation-ray hit locations on the surfaces of a $20 \times 20 \times 20 \, \mathrm{m}^3$ cubic enclosure for varying numbers of random rotations of the angular distribution, with the rays emitted from a spherical source (radius = 1 m) at the enclosure center. Each color represents the rays generated during one random rotation. Increasing the number of rotations improves the spatial uniformity of the surface coverage and reduces directional bias associated with a fixed angular discretization.}
+\label{fig:rad_rand_rot_cube_hit}
+\end{center}
+\end{figure}
+
 
 \subsection{Spatial Discretization in Cut Cells}
 

Manuals/FDS_Technical_Reference_Guide/Radiation_Chapter.tex

@@ -10444,7 +10444,7 @@ \subsection{Flame Height and Flame Tilt}
 The ``flame height,'' and to some extent the ``flame tilt,'' are common quantities used to assess the hazard from a large fire. FDS does not output these quantities directly, but there is a combination of devices and controls that yield these values. This is all best explained with an example. Suppose you are simulating a fire in a domain that is 10~m high with a grid resolution of 0.2~m. The fire is centered about the point $(x,y)=(0,0)$ at $z=0$. The following input lines compute the flame height and tilt:
 \begin{lstlisting}
 &DEVC ID='HRR', Z_ID='z', XBP=0.0,0.0,0.0,0.0,0.1,9.9, QUANTITY='HRRPUV',
-SPATIAL_STATISTIC='VOLUME INTEGRAL', DX=5, DY=5, DZ=0.1, POINTS=100,
+SPATIAL_STATISTIC='VOLUME INTEGRAL', DX=5, DY=5, DZ=0.1, POINTS=50,
 STATISTICS_START=10. /
 &CTRL ID='H', FUNCTION_TYPE='PERCENTILE', INPUT_ID='HRR', PERCENTILE=0.97 /
 &DEVC ID='x_max', XB=-5,5,-5,5,9.0,9.2, QUANTITY='TEMPERATURE',

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -3629,6 +3629,7 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT)
       WRITE(LU_OUTPUT,'(A,ES10.3)') '     Packing Ratio               ', SF%VEG_LSET_BETA
       WRITE(LU_OUTPUT,'(A,ES10.3)') '     Surface Area/Volume (1/m)   ', SF%VEG_LSET_SIGMA*100.  ! Convert from 1/cm to 1/m
       WRITE(LU_OUTPUT,'(A,ES10.3)') '     Fuel Depth (m)              ', SF%VEG_LSET_HT
+      WRITE(LU_OUTPUT,'(A,ES10.3)') '     Fuel Load (kg/m2)           ', SF%VEG_LSET_SURF_LOAD
    ENDIF
 
 ENDDO SURFLOOP

Source/dump.f90

@@ -8326,15 +8326,17 @@ SUBROUTINE READ_SURF(QUICK_READ)
    ENDIF
 
    ! Level set vegetation fire spread specific
+   ! VEG_LSET_FUEL_INDEX: -1=undefined, 0=custom Rothermel from layer props, 1-13=Albini fuel models
 
    VEG_LSET_SPREAD = .FALSE.
-   IF (VEG_LSET_IGNITE_TIME < 1.E6_EB .OR. VEG_LSET_FUEL_INDEX>0 .OR. VEG_LSET_ROS_00>0._EB) VEG_LSET_SPREAD = .TRUE.
+   ! SURF_LOAD without ROS_00/index implies custom fuel (index 0)
+   IF (VEG_LSET_FUEL_INDEX<0 .AND. VEG_LSET_ROS_00<=0._EB .AND. VEG_LSET_SURF_LOAD>0._EB) VEG_LSET_FUEL_INDEX = 0
+   IF (VEG_LSET_IGNITE_TIME < 1.E6_EB .OR. VEG_LSET_FUEL_INDEX>=0 .OR. VEG_LSET_ROS_00>0._EB) VEG_LSET_SPREAD = .TRUE.
    IF (VEG_LSET_SPREAD .AND. LEVEL_SET_MODE==0) THEN
       WRITE(MESSAGE,'(A,A,A)') 'ERROR(305): SURF ',TRIM(ID),' indicates a level set simulation, but LEVEL_SET_MODE not set on MISC.'
       CALL SHUTDOWN(MESSAGE) ; RETURN
    ENDIF
-   IF (VEG_LSET_FUEL_INDEX>0 .AND. LEVEL_SET_COUPLED_FIRE) HRRPUA = 1._EB  ! HRRPUA to be set properly later
-   IF (VEG_LSET_ROS_00    >0 .AND. LEVEL_SET_COUPLED_FIRE) HRRPUA = 1._EB
+   IF ((VEG_LSET_FUEL_INDEX>=0 .OR. VEG_LSET_ROS_00>0._EB) .AND. LEVEL_SET_COUPLED_FIRE) HRRPUA = 1._EB  ! set properly later
 
    SF%VEG_LSET_SPREAD       = VEG_LSET_SPREAD
    SF%VEG_LSET_ROS_00       = VEG_LSET_ROS_00       ! no-wind, no-slope RoS (m/s), Rothermel model
@@ -8376,8 +8378,8 @@ SUBROUTINE READ_SURF(QUICK_READ)
       END SELECT
    ENDIF
 
-   ! Set defaults for custom fuel model (for indexed fuels they are either user-specified or calculated later)
-   IF (SF%VEG_LSET_FUEL_INDEX==0._EB) THEN
+   ! Defaults for undefined fuel properties (set elsewhere for indexed fuels)
+   IF (SF%VEG_LSET_FUEL_INDEX<=0) THEN
       IF (SF%VEG_LSET_SIGMA<0._EB) SF%VEG_LSET_SIGMA = 50._EB
       IF (SF%VEG_LSET_HT<0._EB) SF%VEG_LSET_HT = 0.2_EB
       IF (SF%VEG_LSET_SURF_LOAD<0._EB) SF%VEG_LSET_SURF_LOAD = 1.0_EB
@@ -9306,7 +9308,7 @@ SUBROUTINE SET_SURF_DEFAULTS
 VEG_LSET_M100           = 0.05_EB
 VEG_LSET_MLW            = 0.70_EB
 VEG_LSET_MLH            = 0.70_EB
-VEG_LSET_FUEL_INDEX     = 0
+VEG_LSET_FUEL_INDEX     = -1  ! -1=undefined, 0=custom Rothermel, 1-13=Albini models
 VEG_LSET_SURF_LOAD      = -1.0_EB !kg/m^2
 VEG_LSET_FIREBASE_TIME  = -1.0_EB
 VEG_LSET_CHAR_FRACTION  = 0.20_EB

Source/read.f90

@@ -46,17 +46,18 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_1(NM)
 
 M => MESHES(NM)
 
-! Loop through all SURFace types and find level set cases that need a calculated RoS
+! Loop through all SURFace types and find level set cases that need a calculated RoS.
+! VEG_LSET_FUEL_INDEX: -1=undefined (user ROS_00), 0=custom Rothermel, 1-13=Albini models.
 
 DO SURF_INDEX=0,N_SURF
    SF => SURFACE(SURF_INDEX)
-   IF (SF%VEG_LSET_SPREAD .AND. SF%VEG_LSET_FUEL_INDEX>0) THEN
+   IF (.NOT.SF%VEG_LSET_SPREAD) CYCLE
+   IF (SF%VEG_LSET_FUEL_INDEX>=0) THEN
       SF%VEG_LSET_ROS_00 = ROS_NO_WIND_NO_SLOPE(SF%VEG_LSET_FUEL_INDEX,SURF_INDEX)
-   ENDIF
-   IF (SF%VEG_LSET_SPREAD .AND. SF%VEG_LSET_FUEL_INDEX==0) THEN
-     SF%BURN_DURATION = SF%VEG_LSET_FIREBASE_TIME
-     IF (LEVEL_SET_COUPLED_FIRE) SF%MASS_FLUX(REACTION(1)%FUEL_SMIX_INDEX) = &
-       (1._EB-SF%VEG_LSET_CHAR_FRACTION)*SF%VEG_LSET_SURF_LOAD/SF%VEG_LSET_FIREBASE_TIME
+   ELSE
+      SF%BURN_DURATION = SF%VEG_LSET_FIREBASE_TIME
+      IF (LEVEL_SET_COUPLED_FIRE) SF%MASS_FLUX(REACTION(1)%FUEL_SMIX_INDEX) = &
+        (1._EB-SF%VEG_LSET_CHAR_FRACTION)*SF%VEG_LSET_SURF_LOAD/SF%VEG_LSET_FIREBASE_TIME
    ENDIF
 ENDDO
 
@@ -1010,6 +1011,11 @@ REAL(EB) FUNCTION ROS_NO_WIND_NO_SLOPE(ROTHERMEL_FUEL_INDEX,SURF_INDEX)
 mlh = SF%VEG_LSET_MLH
 
 SELECT CASE(ROTHERMEL_FUEL_INDEX)
+   CASE(0)  ! Custom homogeneous fuel from SURF layer properties
+      w0d1=SF%VEG_LSET_SURF_LOAD ; w0d2=0._EB ; w0d3=0._EB ; w0lh=0._EB ; w0lw=0._EB
+      svd1=SF%VEG_LSET_SIGMA*100._EB  ! stored as 1/cm; Rothermel formulas use 1/m
+      svd2=358._EB ; svd3=98._EB ; svlh=4921._EB ; svlw=4921._EB
+      mx=0.15_EB ; depth=SF%VEG_LSET_HT ; rhop=512._EB ; heat=18607._EB ; st=0.0555_EB ; se=0.01_EB
    CASE(1)  ! 'Short Grass'
       w0d1=0.1659     ; w0d2=0.        ; w0d3=0.        ; w0lh=0.        ; w0lw=0.     ! dry mass per unit area (kg/m2)
       svd1=11483.     ; svd2=358.      ; svd3=98.       ; svlh=4921.     ; svlw=4921.  ! surface area to volume (1/m)
@@ -1076,6 +1082,8 @@ REAL(EB) FUNCTION ROS_NO_WIND_NO_SLOPE(ROTHERMEL_FUEL_INDEX,SURF_INDEX)
    w0lh = SF%VEG_LSET_SURF_LOAD/w0*w0lh
    w0lw = SF%VEG_LSET_SURF_LOAD/w0*w0lw
    w0 = (w0d1 + w0d2 + w0d3 + w0lh + w0lw)
+ELSE
+   SF%VEG_LSET_SURF_LOAD = w0
 ENDIF
 
 ! Auxiliary functions