Merge pull request #15642 from rmcdermo/master

FDS Tech Guide: add section describing FLUX_LIMITER_MW_CORRECTION

Author: rmcdermo

Manuals/FDS_Technical_Reference_Guide/Equation_Chapter.tex

@@ -149,9 +149,6 @@ \section{Mass and Species Transport}
 \be \dod{\rho}{t} + \nabla\!\cdot (\rho \bu)  =  \dm_{\rm b}'''  \label{mass} \ee
 because $\sum Z_\alpha=1$ and $\sum \dm_\alpha''' = 0$ and $\sum \dm_{\rm b,\alpha}'''=\dm_{\rm b}'''$, by definition, and because it is assumed that $\sum \rho D_\alpha \nabla Z_\alpha = 0$. This last assertion is not true, in general.  The diffusive flux for the most abundant local species is corrected to enforce the constraint.
 
-
-
-
 \paragraph{Enforcing Realizability} ~\\
 
 \noindent Realizability of species mass fractions requires $Y_\alpha \ge 0$ for all $\alpha$ and $\sum Y_\alpha = 1$.  Note that this is more restrictive than the boundedness constraint, which simply requires $0 \le Y_\alpha \le 1$.
@@ -160,6 +157,10 @@ \section{Mass and Species Transport}
 
 With this approach we must take care to ensure $\sum \rho D_\alpha \nabla Z_\alpha = 0$.  Our strategy is to absorb any errors in diffusive transport into the most abundant species \emph{locally}.  That is, for a given cell face we set $\rho D_m \nabla Z_m = -\sum_{\alpha\ne m} \rho D_\alpha \nabla Z_\alpha$, where $m$ is the most abundant species adjacent to that face.  Note that since FDS is typically used as an LES code mass transport by molecular diffusion may be two or three orders of magnitude less than turbulent transport, which uses the same turbulent diffusion coefficient for all species. Therefore, the errors in summation of the diffusive fluxes tend to be small.
 
+\paragraph{Enforcing the Equation of State} ~\\
+
+\noindent Given the strategy of solving $N_s$ transport equations and obtaining $\rho = \sum_{\alpha=1}^{N_s} (\rho Y)_\alpha$ to enforce realizability, we must take care to ensure that we maintain the equation of state at cell faces after flux limited interpolation.  In Sec.~\ref{sec_flux_limiters}, we describe a correction scheme that enforces the equation of state at cell faces and thus maintains isothermal flows for multi-component mixtures with variable molecular weights.
+
 \section{Low Mach Number Approximation}
 
 For low speed applications like fire, Rehm and Baum~\cite{Rehm:1} observed that the spatially and temporally resolved pressure, $p$, can be decomposed into a ``background'' pressure, $\bp(z,t)$, plus a perturbation, $\tp(x,y,z,t)$, with only the background pressure retained in the equation of state (ideal gas law):

Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex

@@ -78,7 +78,7 @@ \subsection{Flux Limiters}
 
 For uniform flow velocity, a fundamental property of the exact solution to the equations governing scalar transport is that the total variation of the scalar field (the sum of the absolute values of the scalar differences between neighboring cells) is either preserved or diminished (never increased).  In other words, no new extrema are created.  Numerical schemes which preserve this property are referred to as total variation diminishing (TVD) schemes.  The practical importance of using a TVD scheme for fire modeling is that such a scheme is able to accurately track coherent vortex structure in turbulent flames and does not develop spurious reaction zones.
 
-FDS employs two second-order TVD schemes as options for scalar transport: Superbee and CHARM.  Superbee \cite{Roe:1986} is recommended for LES because it more accurately preserves the scalar variance for coarse grid solutions that are not expected to be smooth.  Due to the gradient steepening applied in Superbee, however, the convergence degrades at small grid spacing for smooth solutions (the method will revert to a stair-step pattern instead of the exact solution).  CHARM \cite{Zhou:1995}, though slightly more dissipative than Superbee, is convergent, and is therefore the better choice for DNS calculations where the flame front is well resolved.
+FDS employs two second-order TVD schemes as options for scalar transport: Superbee and CHARM.  Superbee \cite{Roe:1986} is recommended for VLES because it more accurately preserves the scalar variance for coarse grid solutions that are not expected to be smooth.  Due to the gradient steepening applied in Superbee, however, the convergence degrades at small grid spacing for smooth solutions (the method will revert to a stair-step pattern instead of the exact solution).  CHARM \cite{Zhou:1995}, though slightly more dissipative than Superbee, is convergent, and is therefore the better choice for LES and DNS calculations.
 
 To illustrate how flux limiters are applied to the scalar transport equations, below we discretize the advection terms in Eq.~(\ref{species}) in one dimension:
 \be  \frac{(\rho Z)_{i}^* - (\rho Z)_{i}^n}{\dt}
@@ -94,12 +94,12 @@ \subsection{Flux Limiters}
 \begin{table}[H]
 \begin{center}
 \begin{tabular}{lc}
-Flux Limiter                           & $B(r)$                         \\
+Flux Limiter                                    & $B(r)$                         \\
 \hline
-Central Difference                     & 1                              \\
-Godunov                                & 0                              \\
-Superbee \cite{Roe:1986} (LES default) & $\max(0,\min(2r,1),\min(r,2))$ \\
-CHARM \cite{Zhou:1995} (DNS default)   & $s(3s+1)/(s+1)^2$; $s=1/r$     \\
+Central Difference                              & 1                              \\
+Godunov                                         & 0                              \\
+Superbee \cite{Roe:1986} (VLES default)         & $\max(0,\min(2r,1),\min(r,2))$ \\
+CHARM \cite{Zhou:1995} (LES and DNS default)    & $s(3s+1)/(s+1)^2$; $s=1/r$     \\
 \end{tabular}
 \end{center}
 \end{table}
@@ -123,21 +123,38 @@ \subsubsection{Notes on Implementation}
 
 The Central Difference and Godunov limiters are included for completeness, debugging, and educational purposes.  These schemes have little utility for typical FDS applications.
 
-% \subsubsection{Dealing with Variable Molecular Weights}
+\subsubsection{Dealing with Variable Molecular Weights}
 
-%% leave this here for a moment as a reminder to write up the constant limiter coefficient method we are now using.
+Handling multi-component mixtures with variable molecular weights requires special care when constructing fluxes, otherwise we may violate the equation of state.  The problem is easiest to illustrate by considering an isothermal (not necessarily constant density) flow.  Maintaining isothermal flow requires
+\begin{equation}
+T = \frac{\overline{W} \bar{p}}{\rho R} = \frac{\bar{p}}{R \rho \sum_\alpha \frac{Z_\alpha}{W_\alpha}} = \frac{\bar{p}}{R \sum_\alpha \frac{\rho Z_\alpha}{W_\alpha}}
+\end{equation}
+to be constant and uniform at all cells and faces. Therefore, with $\bar{p}$ and $R$ constant and uniform, we must maintain
+\begin{equation}
+\label{eq:rho_mw}
+\sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}
+\end{equation}
+The above condition is automatically satisfied in the cases of using Godunov or Central differencing or in the case of binary flow (two species).  However, if we apply a second-order flux limiter, such as Superbee or CHARM, independently to each species in a multi-component (three or more species) flow with variable molecular weights, then this condition is easily violated at a cell face.
+
+For example, consider the case of a multi-component, constant density flow.  A flux limited face value for species $\alpha$ is a sum of the values in adjacent cells times a coefficient $c_{i,\alpha}$ for each cell $i$.  At a face we require have $\sum_\alpha \sum_i c_{i,\alpha} (\rho Z_\alpha)_i = \rho = \mbox{constant}$.  The only way to guarantee $\rho$ is constant is if the $c_{i,\alpha}$ are the same for each $\alpha$, that is, $c_{i,\alpha} = c_i$.  Rearranging then gives $\sum_i c_i \sum_\alpha (\rho Z_\alpha)_i = \sum_i c_i \rho_i = \rho$, which just says that we need to obey the constraint $\sum_i c_i = 1$, which is true by construction of the limiter.
 
-% Maintaining isothermal flow requires
-% \begin{equation}
-% T = \frac{\overline{W} \bar{p}}{\rho R} = \frac{\bar{p}}{R \rho \sum_\alpha \frac{Z_\alpha}{W_\alpha}} = \frac{\bar{p}}{R \sum_\alpha \frac{\rho Z_\alpha}{W_\alpha}}
-% \end{equation}
-% to be constant and uniform at all cells and faces. Therefore, with $\bar{p}$ and $R$ constant and uniform, we must maintain
-% \begin{equation}
-% \label{eq:rho_mw}
-% \sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}
-% \end{equation}
+As mentioned above, constant limiter coefficients are only possible with simple schemes like Godunov (low order) or central (not TVD).  In practice, applying a second-order limiter to all species independently and then taking the worst case also reduces to a low order scheme.  Therefore, it is important to maintain the independent second-order limiters for the individual species.  An algorithm to achieve this goal and satisfy the equation of state is described next.
 
-% The above condition is automatically satisfied in the cases of using Godunov or Central differencing or in the case of binary flow (two species).  However, if we apply a second-order flux limiter, such as Superbee or CHARM, independently to each species in a multi-component (three or more species) flow with variable molecular weights, then this condition is easily violated.
+\subsubsection{Molecular Weight Correction Algorithm}
+
+Whether the flow is isothermal or not, satisfying Eq.~(\ref{eq:rho_mw}) is required to obey the equation of state.  The scheme introduced here uses a flux limiter to find a face value of $\rho/\overline{W}$ and then corrects the most abundant species face density to satisfy Eq.~(\ref{eq:rho_mw}).
+
+The following algorithm computes flux-limited species face densities for multi-component mixtures with variable molecular weights and preserves the equation of state, Eq.~(\ref{eq:rho_mw}), at cell faces.  For simplicity, the following considers a 1-D domain with cell face between cells $i$ and $i+1$.
+\begin{enumerate}
+\item Compute flux limited face values as usual, $\overline{(\rho Z_\alpha)}^{FL}_{i+1/2}$, for $\alpha=1..N_s$.
+\item Also, for each face, determine the {\tt maxloc} of the flux-limited $\overline{(\rho Z_\alpha)}^{FL}_{i+1/2}$ (the most abundant species on the face) and call it $\gamma$.  We will absorb the error into this species.
+\item For each cell, compute $\displaystyle \frac{\rho}{\overline{W}}$ from Eq.~(\ref{eq:rho_mw}).
+\item For each cell face, compute $\displaystyle \overline{\left(\frac{\rho}{\overline{W}}\right)}^{FL}_{i+1/2}$.
+\item Last, determine the face component density for $\gamma$ as
+\begin{equation}
+\overline{(\rho Z_\gamma)}^{FL}_{i+1/2} = W_\gamma \left[\overline{\left(\frac{\rho}{\overline{W}}\right)}^{FL}_{i+1/2} - \sum_{\alpha\ne\gamma} \frac{\overline{(\rho Z_{\alpha})}^{FL}_{i+1/2}}{W_\alpha} \right]
+\end{equation}
+\end{enumerate}
 
 
 \subsection{Time Splitting for Mass Source Terms}