\documentstyle[11pt]{article}
\begin{document}

\begin{center}
{\bf The Three-Dimensional Compressible Navier-Stokes Equations}
\end{center}

\noindent The Navier-Stokes equations in conservative form are given by

\begin{eqnarray}
\partial_t^*\,\rho^* + \partial_j^* \left( \rho^*\,u_j^* \right) & = & 0\ , 
\nonumber \\
\partial_t^* \left( \rho^*\,u_i^* \right) + \partial_j^* \left\lbrack \rho^* u_i^* 
u_j^* + p^* \delta_{ij} - \tau^*_{ij} \right\rbrack & = & \rho^* F_i^*\ , 
\nonumber \\
\partial_t^* (\rho^* E^*) + \partial_j^* \left\lbrack \rho^* u_j^* H^* - u_i^* 
\tau^*_{ij} - k^* \partial_j^* T^* \right\rbrack & = & \rho^* u_i^* F_i^* + 
\dot{q}^*\ ,
\end{eqnarray}

\noindent where the superscript `*' denotes dimensional quantities and

\begin{equation}
\tau^*_{ij} = \mu^* \left\lbrack \partial_j^* u_i^* + \partial_i^* u_j^* 
\right\rbrack + \lambda^*\,\delta_{ij}\,\partial^*_k\,u^*_k\ ,
\end{equation}

\begin{description}
\item{$\rho^*$:} density
\item{$u_i^*$:} $i^{th}$ velocity component, $i=1,3$ 
\item{$p^*$:} pressure
\item{$F_i^*$:} $i^{th}$ component of the external force
\item{$E^*$:} total internal energy $ = e^* + \frac{u_i^*\,u_i^*}{2}$, where $e^*$ is
the specific internal energy
\item{$H^*$:} total enthalpy $ = h^* + \frac{u_i^*\,u_i^*}{2}$, where $h^*$ is the
specific enthalpy
\item{$T^*$:} temperature
\item{$k^*$:} thermal conductivity
\item{$\mu^*$:} absolute (or dynamic) viscosity
\item{$\lambda^*$:} second coefficient of viscosity ($= - \frac{2}{3}\mu^*$ upon
invoking Stokes assumption)
\item{$\dot{q}^*$:} volumetric heat generation
\end{description}

The specific enthalpy is related to the specific internal energy by

\begin{equation}
h^* = e^* + \frac{p^*}{\rho^*} \Rightarrow H^* = E^* + \frac{p^*}{\rho^*}\ .
\end{equation}

A perfect gas is assumed, in which case,

\begin{equation}
e^* = E^* - \frac{u_i^*\,u_i^*}{2} = C_v^* T^* \Rightarrow T^* = \frac{1}{C_v^*}
\left\lbrack E^* - \frac{u_i^*\,u_i^*}{2} \right\rbrack\ .
\end{equation}
\vspace{10pt}

\noindent Furthermore, for a perfect gas, $p^* = \rho^* R^* T^*$. Hence,

\begin{equation}
p^* = \rho^* \frac{R^*}{C_v^*} \left\lbrack E^* - \frac{u_i^*\,u_i^*}{2}
\right\rbrack = (\gamma -1) \rho^* \left\lbrack E^* - \frac{u_i^*\,u_i^*}{2} 
\right\rbrack\ ,
\end{equation}

\noindent where

\begin{description}
\item{$C_v^*$:} specific heat at constant volume
\item{$R^*$:} gas constant
\item{$\gamma$:} ratio of specific heats $= \frac{C_p^*}{C_v^*}$
\item{$C_p^*$:} specific heat at constant pressure
\end{description}
\vspace{15pt}

\noindent{\bf Non-dimensional Form of the Equations}

We non-dimensionalize the equations by defining the following reference
values:

\begin{description}
\item{$U_\infty^*$:} reference velocity
\item{$L_\infty^*$:} reference length scale
\item{$\rho_\infty^*$:} reference density
\item{$\mu_\infty^*$:} reference absolute viscosity
\item{$k_\infty^*$:} reference thermal conductivity
\item{$C_{v\infty}^*$:} reference specific heat
\end{description}
\vspace{10pt}

All non-dimensional quantities are denoted without the superscript `*'.
Hence,

\begin{eqnarray}
u_i & = & \frac{u_i^*}{u_\infty^*} \ ;\hspace{1em}
x_i = \frac{x_i^*}{L_\infty^*} \ ;\hspace{1em}
\rho = \frac{\rho^*}{\rho_\infty^*}= ; \nonumber \\
\mu & = & \frac{\mu^*}{\mu_\infty^*} \ ;\hspace{1em} \lambda =
\frac{\lambda^*}{\mu_\infty^*}\ ;\hspace{1em} k = \frac{k^*}{k_\infty^*}\
;\hspace{1em} C_v = \frac{C_v^*}{C_{v\infty}^*}\ .
\end{eqnarray}

\noindent In addition, we non-dimensionalize the other quantities in
equations~(1)--(4) as follows:

\begin{eqnarray}
t & = & \frac{t^*U_\infty^*}{L^*} \ ;\hspace{1em}
p = \frac{p^*}{\rho_\infty^*\,U_\infty^{*2}} \ ;\hspace{1em}
F_i = \frac{F_i^*}{U_\infty^{*2}/L_\infty^*} \ ;\hspace{1em}
E = \frac{E^*}{U_\infty^{*2}} \nonumber \\
H & = & \frac{H^*}{U_\infty^{*2}} \ ;\hspace{1em}
\dot{q} = \frac{\dot{q}^*}{\rho_\infty^*\, U_\infty^{*3}\,/\, L_\infty^*}
\ ;\hspace{1em}
\tau_{ij} = \frac{\tau_{ij}^*}{\left( U_\infty^*\, \mu_\infty^*
\right)\,/\, L_\infty^*}\ .
\end{eqnarray}

Upon substituting equations~(6) and (7) into equations~(1)--(2) and
using equation~(4), it may be easily shown that

\begin{eqnarray}
\partial_t \rho + \partial_j \left( \rho u_j \right) & = & 0\ , \nonumber \\
&& \nonumber \\
\partial_t \left( \rho u_i \right) + \partial_j \left\lbrack \rho u_i u_j + p
\delta_{ij} - \frac{1}{\rm Re}\, \tau_{ij} \right\rbrack & = & \rho F_i\ , \nonumber \\
&& \nonumber \\
\partial_t (\rho E) + \partial_j \bigg[ \rho u_j H - \frac{1}{\rm Re}\, u_i
\tau_{ij} & - & \nonumber \\
\frac{\gamma}{\rm Re\,Pr}\, k \partial_j \left\lbrace\frac{1}{C_v} 
\left( E - \frac{u_i u_i}{2} \right) \right\rbrace \bigg] & = & \rho u_i F_i + 
\dot{q}\ ,
\end{eqnarray}

\noindent where
 
\begin{equation}
\tau_{ij} = \mu \left( \partial_j u_i + \partial_i u_j \right) -
\frac{2}{3} \mu \delta_{ij} \partial_k u_k\ ,
\end{equation}

\noindent and

\[
\mbox{Re\ }\ = \frac{\rho_\infty^* u_\infty^*
L_\infty^*}{\mu_\infty^*}\ , 
\]
\[
\mbox{Pr\ }\ = \frac{\mu_\infty^* C_{p\infty}^*}{k_\infty^*}
= \frac{\mu_\infty^* C_{v\infty}^* \gamma}{k_\infty^*}\ .
\]

\noindent Here, Re and Pr denote the Reynolds and Prandtl numbers,
respectively.  Finally, the non-dimensional form of equations~(3) and
(4) are given by

\begin{eqnarray}
H & = & E + \frac{p}{\rho}\ , \\
p & = & (\gamma - 1) \rho \left\lbrack E - \frac{u_i u_i}{2} \right\rbrack\
.
\end{eqnarray}

\noindent{\bf Note:}\ The Stokes' assumption has been invoked in equation~(9).
\vspace{15pt}

\begin{center}
{\bf Three-Dimensional Compressible Navier-Stokes Equations and the NAS
Benchmark}
\end{center}

\noindent Let the system of equations~(8) be written as

\begin{equation}
\frac{\partial {\bf q}}{\partial t} +
\frac{\partial {\bf E}}{\partial x_1} +
\frac{\partial {\bf F}}{\partial x_2} +
\frac{\partial {\bf G}}{\partial x_3} =
\frac{\partial {\bf U}}{\partial x_1} +
\frac{\partial {\bf V}}{\partial x_2} +
\frac{\partial {\bf W}}{\partial x_3} +  {\bf H}\ ,
\end{equation}

\noindent where

\[
{\bf q} = \left\lbrack \begin{array}{l} \rho \\
\rho u_1 \\
\rho u_2 \\
\rho u_3 \\
\rho E \\ \end{array}
\right\rbrack ; 
\]
\[
{\bf E} = \left\lbrack \begin{array}{l} \rho u_1 \\
\rho u_1^2 + p \\
\rho u_1 u_2 \\
\rho u_1 u_3 \\
u_1 (\rho E + p) \\ \end{array}
\right\rbrack
; \qquad {\bf F} = \left\lbrack \begin{array}{l} \rho u_2 \\
\rho u_1 u_2 \\
\rho u_2^2 + p \\
\rho u_2 u_3 \\
u_2 (\rho E + p) \\ \end{array}
\right\rbrack ;
\qquad {\bf G} = \left\lbrack \begin{array}{l} \rho u_3 \\
\rho u_1 u_3 \\
\rho u_2 u_3 \\
\rho u_3^2 + p \\
u_3 (\rho E + p) \\ \end{array}
\right\rbrack ;
\]

\[
{\bf U} = {\rm Re}^{-1} \left\lbrack \begin{array}{l} 0 \\
\tau_{11} \\
\tau_{21} \\
\tau_{31} \\
\alpha_1 \\ \end{array}
\right\rbrack
; \qquad {\bf V} = {\rm Re}^{-1} \left\lbrack \begin{array}{l} 0 \\
\tau_{12} \\
\tau_{22} \\
\tau_{32} \\
\alpha_2 \\ \end{array}
\right\rbrack
; \qquad {\bf W} = {\rm Re}^{-1} \left\lbrack \begin{array}{l} 0 \\
\tau_{13} \\
\tau_{23} \\
\tau_{33} \\
\alpha_3 \\ \end{array}
\right\rbrack ;
\]

\noindent and

\[
{\bf H} = \left\lbrack \begin{array}{l} 0 \\
\rho F_1 \\
\rho F_2 \\
\rho F_3 \\
\rho u_k F_k + \dot{q} \\ \end{array}
\right\rbrack .
\]

\noindent The quantities $\alpha_i$, $i=1$, 3 are given by

\[
\alpha_1 = u_k \tau_{k1} + \frac{\gamma k}{\rm Pr} \partial_1
\left\lbrace \frac{1}{C_v} \left( E - \frac{u_i u_i}{2}\right) \right\rbrace\ , 
\]
\[
\alpha_2 = u_k \tau_{k2} + \frac{\gamma k}{\rm Pr} \partial_2
\left\lbrace \frac{1}{C_v} \left( E - \frac{u_i u_i}{2}\right) \right\rbrace\ , 
\]
\[
\alpha_3 = u_k \tau_{k3} + \frac{\gamma k}{\rm Pr} \partial_3
\left\lbrace \frac{1}{C_v} \left( E - \frac{u_i u_i}{2}\right) \right\rbrace\ . 
\]
\vspace{10pt}

Let
\[ {\bf U} = {\bf U}_I + {\bf U}_E\ , \]
\[ {\bf V} = {\bf V}_I + {\bf V}_E\ , \]
\[ {\bf W} = {\bf W}_I + {\bf W}_E\ , \]

\noindent so that
$\frac{\partial {\bf U}_I}{\partial x_1}$\ , 
$\frac{\partial {\bf V}_I}{\partial x_2}$\ ,  and
$\frac{\partial {\bf W}_I}{\partial x_3}$ involve no cross derivatives, and
$\frac{\partial {\bf U}_E}{\partial x_1}$\ , 
$\frac{\partial {\bf V}_E}{\partial x_2}$\ , 
$\frac{\partial {\bf W}_E}{\partial x_3}$ do involve cross derivatives.
\vspace{10pt}

As demonstrated in the development of the Beam and Warming algorithm, only
those terms which do not involve cross-derivatives are handled implicitly.
Cross-derivative terms are handled explicitly.  The focus of the NAS
benchmark is on the implicit solution algorithm, and not on the evaluation of
the right-hand side.  Therefore, in effect, a solution of the following
equation is sought:

\begin{equation}
\frac{\partial {\bf q}}{\partial t} +
\frac{\partial {\bf E}}{\partial x_1} +
\frac{\partial {\bf F}}{\partial x_2} +
\frac{\partial {\bf G}}{\partial x_3} -
\frac{\partial {\bf U}_I}{\partial x_1} -
\frac{\partial {\bf V}_I}{\partial x_2} -
\frac{\partial {\bf W}_I}{\partial x_3} =  {\bf H}^*\ ,
\end{equation}

\noindent where

\begin{equation}
{\bf H}^* = {\bf H} + 
\frac{\partial {\bf U}_E}{\partial x_1} +
\frac{\partial {\bf V}_E}{\partial x_2} +
\frac{\partial {\bf W}_E}{\partial x_3}\ . 
\end{equation}

In a real CFD application, equation~(4) would be evaluated at each time
step (for reference, see equation~(22) of BW).  However, in the NAS
benchmark, the individual contributions to ${\bf H}^*$ are ignored, but
in order to guarantee the existence of a steady-state solution, the
function ${\bf H}^*$ is specified in the following manner.  Let ${\bf
q}^*$ denote the steady-state solution (which, for the NAS benchmark,
is given).  Then at steady state,

\begin{equation}
{\bf H}^* = \frac{\partial {\bf E}^*}{\partial x_1} +
\frac{\partial {\bf F}^*}{\partial x_2} +
\frac{\partial {\bf G}^*}{\partial x_3} -
\frac{\partial {\bf U}_I^*}{\partial x_1} -
\frac{\partial {\bf V}_I^*}{\partial x_2} -
\frac{\partial {\bf W}_I^*}{\partial x_3}\ . 
\end{equation}

\noindent If ${\bf H}^*$ is specified according to the above identity,
a steady-state solution of equation~(13) exists.  Thus, for the NAS
benchmark, we need to consider only those terms that do not involve
cross-derivatives.
\vspace{10pt}

First consider the solution ${\bf q}$ and the inviscid fluxes (for
reference, see page~46 of NAS benchmark).
\vspace{10pt}

\begin{tabular}{lcl}
{\bf Present Notation}&&{\bf NAS Benchmark Notation} \\
&& \\
${\bf q} = \left\lbrack \begin{array}{l} \rho \\
\rho u_1 \\
\rho u_2 \\
\rho u_3 \\
\rho E \\ \end{array} \right\rbrack$
& \hspace{1em} $\equiv$ \hspace{1em} &
${\bf\phantom{-}U} = \left\lbrack \begin{array}{l} u^{(1)} \\
u^{(2)} \\
u^{(3)} \\
u^{(4)} \\
u^{(5)} \\ \end{array} \right\rbrack$ \\
&& \\
${\bf E} = \left\lbrack \begin{array}{l} \rho u_1 \\
\rho u_1^2 + p \\
\rho u_1 u_2 \\
\rho u_1 u_3 \\
u_1 (\rho E + p) \\ \end{array} \right\rbrack$
& \hspace{1em} $\equiv$ \hspace{1em} &
${\bf -E} = \left\lbrack \begin{array}{l} u^{(2)} \\
\left\lbrack u^{(2)}\right\rbrack^2\,/\,\left\lbrack u^{(1)}\right\rbrack + \phi \\
u^{(2)} u^{(3)}\,/\, u^{(1)} \\
u^{(2)} u^{(4)}\,/\, u^{(1)} \\
u^{(2)}\,/\, u^{(1)} \left\lbrack u^{(5)} + \phi \right\rbrack \\
\end{array} \right\rbrack$ \\
&& \\
${\bf F} = \left\lbrack \begin{array}{l} \rho u_2 \\
\rho u_1 u_2 \\
\rho u_2^2 + p \\
\rho u_2 u_3 \\
u_2 (\rho E + p) \\ \end{array} \right\rbrack$
& \hspace{1em} $\equiv$ \hspace{1em} &
${\bf -F} = \left\lbrack \begin{array}{l} u^{(3)} \\
u^{(2)} u^{(3)}\,/\, u^{(1)} \\
\left\lbrack u^{(3)}\right\rbrack^2\,/\,\left\lbrack u^{(1)}\right\rbrack + \phi \\
u^{(3)} u^{(4)}\,/\, u^{(1)} \\
u^{(3)}\,/\, u^{(1)} \left\lbrack u^{(5)} + \phi \right\rbrack \\
\end{array} \right\rbrack$ \\
&& \\
${\bf G} = \left\lbrack \begin{array}{l} \rho u_3 \\
\rho u_1 u_3 \\
\rho u_2 u_3 \\
\rho u_3^2 + p \\
u_3 (\rho E + p) \\ \end{array} \right\rbrack$
& \hspace{1em} $\equiv$ \hspace{1em} &
${\bf -G} = \left\lbrack \begin{array}{l} u^{(4)} \\
u^{(2)} u^{(4)}\,/\, u^{(1)} \\
u^{(3)} u^{(4)}\,/\, u^{(1)} \\
\left\lbrack u^{(4)}\right\rbrack^2\,/\,\left\lbrack u^{(1)}\right\rbrack + \phi \\
u^{(4)}\,/\, u^{(1)} \left\lbrack u^{(5)} + \phi \right\rbrack \\
\end{array} \right\rbrack$ \\
\end{tabular}

\[
p = (\gamma -1) \left\lbrack \rho E - \frac{\rho u_i u_i}{2} \right\rbrack \equiv
\]
\[
\phi = k_2 \left\lbrack u^{(5)} - \frac{1}{2} \left\lbrace 
\frac{\left\lbrack u^{(2)}\right\rbrack^2 + \left\lbrack u^{(3)}\right\rbrack^2 +
\left\lbrack u^{(4)}\right\rbrack^2}{u^{(1)}} \right\rbrace \right\rbrack\ ,
\]
\vspace{10pt}

\begin{equation}
\Rightarrow k_2 = \gamma -1\ .
\end{equation}

Next, consider the momentum equations. It follows from equation~(9) that 

\begin{equation}
\partial_j\,\tau_{ij} = \left\lbrack \begin{array}{l}
\partial_1 \left( \frac{4}{3} \mu \partial_1 u_1 \right) + \partial_2 
\left( \mu \partial_2 u_1 \right) \phantom{\frac{4}{3}} + \partial_3 \left( \mu 
\partial_3 u_1 \right) \phantom{\frac{4}{3}} + \mbox{c.d.} \\
\partial_1 \left( \mu \partial_1 u_2 \right) \phantom{\frac{4}{3}} + \partial_2 
\left( \frac{4}{3} \mu \partial_2 u_2 \right) + \partial_3 \left( \mu 
\partial_3 u_2 \right) \phantom{\frac{4}{3}} + \mbox{c.d.} \\
\partial_1 \left( \mu \partial_1 u_3 \right) \phantom{\frac{4}{3}} + \partial_2 
\left( \mu \partial_2 u_3 \right) \phantom{\frac{4}{3}} + \partial_3 \left( 
\frac{4}{3} \mu \partial_3 u_3 \right) + \mbox{c.d.} \\
\end{array} \right\rbrack\ , 
\end{equation}

\noindent where `c.d.' denotes cross-derivative terms.  Next, consider
the energy equation.  Upon substituting for $\tau_{ij}$ from equation~(9), it
may be shown that

\begin{eqnarray}
&& \partial_j \left\lbrack u_i \tau_{ij} + \frac{\gamma k}{\rm Pr}\, \partial_j 
\left\lbrace \frac{1}{C_v} \left( E - \frac{u_i u_i}{2}\right) \right\rbrace 
\right\rbrack \nonumber \\
&= & \partial_1 \left\lbrack \frac{\mu}{2} \left( 1 - \frac{\gamma k}{{\rm Pr} 
\mu C_v} \right) \partial_1 \left( u_1^2 + u_2^2 + u_3^2 \right) + \frac{\mu}{6} 
\partial_1 \left( u_1^2 \right) + \frac{\gamma k}{{\rm Pr} C_v} 
\partial_1 E \right\rbrack \nonumber \\
&+ & \partial_2 \left\lbrack \frac{\mu}{2} \left( 1 - \frac{\gamma k}{{\rm Pr} 
\mu C_v} \right) \partial_2 \left( u_1^2 + u_2^2 + u_3^2 \right) + \frac{\mu}{6} 
\partial_2 \left( u_2^2 \right) + \frac{\gamma k}{{\rm Pr} C_v} 
\partial_2 E \right\rbrack  \nonumber \\
&+ & \partial_3 \left\lbrack \frac{\mu}{2} \left( 1 - \frac{\gamma k}{{\rm Pr} 
\mu C_v} \right) \partial_3 \left( u_1^2 + u_2^2 + u_3^2 \right) + \frac{\mu}{6}
\partial_3 \left( u_3^2 \right) + \frac{\gamma k}{{\rm Pr} C_v} 
\partial_3 E \right\rbrack  + \mbox{c.d.}\ , \nonumber \\
&& \mbox{where}\ C_v\ \mbox{is assumed const.}
\end{eqnarray}
\vspace{10pt}

Now we can construct the functions ${\bf U}_I$, ${\bf V}_I$, and
${\bf W}_I$ using equations~(17) and (18).
\vspace{10pt}

\begin{tabular}{lcl}
{\bf Present Notation}&&{\bf NAS Benchmark Notation} \\
&& \\
${\bf U}_I = {\rm Re}^{-1} = \left\lbrack \begin{array}{l} 0 \\
\\
\frac{4}{3} \mu \partial_1 u_1 \\
\\
\mu \partial_1 u_2 \\
\\
\mu \partial_1 u_3 \\
\\
u_I^5 \\ \end{array} \right\rbrack$
& \hspace{1em} $\equiv$ \hspace{1em} &
$T = \left\lbrack \begin{array}{l} 0 \\
\\
\frac{4}{3} k_3 k_4 \frac{\partial}{\partial\xi} \left( \frac{u^{(2)}}{u^{(1)}}
\right) \\
\\
k_3 k_4 \frac{\partial}{\partial\xi} \left( \frac{u^{(3)}}{u^{(1)}} \right) \\
\\
k_3 k_4 \frac{\partial}{\partial\xi} \left( \frac{u^{(4)}}{u^{(1)}} \right) \\
\\
t^{(5)} \\ \end{array} \right\rbrack$ \\
\end{tabular}
where
\[
u_I^5 = \frac{\mu}{2} \left( 1 - \frac{\gamma k}{{\rm Pr} \mu C_v} \right)
\partial_1 \left( u_1^2 + u_2^2 + u_3^2 \right) + \frac{\mu}{6} \partial_1
\left( u_1^2 \right) + \frac{\gamma k}{{\rm Pr} C_v} \partial_1 E\ ,
\]

\noindent and 

\begin{eqnarray*}
t^{(5)} & = & \frac{1}{2} \left( 1 - k_1 k_5 \right) \frac{\partial}{\partial\xi}
\left\lbrack \frac{\left\lbrack u^{(2)} \right\rbrack^2 + \left\lbrack
u^{(3)} \right\rbrack^2 + \left\lbrack u^{(4)} \right\rbrack^2}{\left\lbrack
u^{(1)} \right\rbrack^2} \right\rbrack \\
& + & \frac{1}{6} \frac{\partial}{\partial\xi}
\left\lbrack \frac{\left\lbrack u^{(2)} \right\rbrack^2}{\left\lbrack u^{(1)}
\right\rbrack^2} \right\rbrack + k_1 k_5 \frac{\partial}{\partial\xi}
\left\lbrack \frac{u^{(5)}}{u^{(1)}} \right\rbrack\ .
\end{eqnarray*}

Similarly, the correspondence between ${\bf V}_I$ and ${\bf V}$, and,
${\bf W}_I$ and ${\bf W}$, may be constructed.
\vspace{10pt}

\noindent{\bf Discrepancies}
\vspace{10pt}

\noindent{\bf (A)\ \ Physical Properties}

\setcounter{equation}{18}

\begin{itemize}
\item Comparison of elements 2--4 (i.e., the momentum equations) indicates
that
\begin{equation}
k_3 k_4 = {\rm Re}^{-1} \mu\ .
\end{equation}
\item Comparison of the energy equation indicates inconsistencies.  The two
forms may be reconciled if $t^{(5)}$ is correctly rewritten as
\begin{eqnarray*}
t^{(5)} & = & k_3 k_4 \bigg\{ \frac{1}{2} \left( 1- k_1 k_5 \right) 
\frac{\partial}{\partial\xi}
\left\lbrack \frac{\left\lbrack u^{(2)} \right\rbrack^2 + \left\lbrack
u^{(3)} \right\rbrack^2 + \left\lbrack u^{(4)} \right\rbrack^2}{\left\lbrack
u^{(1)} \right\rbrack^2} \right\rbrack \\
& + & \frac{1}{6} \frac{\partial}{\partial\xi}
\left\lbrack \frac{\left\lbrack u^{(2)} \right\rbrack^2}{\left\lbrack u^{(1)}
\right\rbrack^2} \right\rbrack + k_1 k_5 \frac{\partial}{\partial\xi}
\left\lbrack \frac{u^{(5)}}{u^{(1)}} \right\rbrack \bigg\}\ ,
\end{eqnarray*}

\setcounter{equation}{19}

where
\begin{equation}
k_1 k_5 = \frac{\gamma k}{{\rm Pr} \mu C_v}\ .
\end{equation}
\end{itemize}

{\bf Note:\ }Although the above factor of $k_3 k_4$ is missing in the
definitions of $t^{(5)}$, $v^{(5)}$, and $w^{(5)}$ (c.f. see pages 47--48
NAS-BM), it is correctly taken into account subsequently in the analytical
evaluation of the various Jacobian matrices (c.f., see $n_{51}\ldots n_{55}$;
$q_{51}\ldots q_{55}$; $s_{51}\ldots s_{55}$ on pages 55--56 NAS-BM).
\vspace{10pt}

\noindent{\bf Discussion of the values of $k_i$, $i=1,5$}
\vspace{10pt}

Although the values of $k_i$ are immaterial insofar as the parallel
implementation of the computational algorithm is concerned, we discuss here
how well these values correspond to those encountered in a ``real'' CFD
application code.

The suggested values in the NAS benchmark are as follows:
\[
k_1 = 1.4 ,\hspace{0.5em}
k_2 = 0.4 ,\hspace{0.5em}
k_3 = 0.1 ,\hspace{0.5em}
k_4 = 1.0 ,\hspace{0.5em}
k_5 = 1.4 .
\]

\noindent However, the above appear in only three groups, viz., $k_{2},
\left( k_1 k_5 \right)$ and $\left( k_3 k_4 \right)$.  The first
constant $k_2$ was found to be $k_2 = \gamma -1$ (c.f. equation~(16)).
Since the ratio of specific heats, $\gamma$, for a diatomic gas is 1.4
(such as air which is composed mainly of $\rm N_2 + O_2$), the above
value for $k_2$ is appropriate.

Next, consider $k_1 k_5$ which is given by equation~(20).
Recall that ${\rm Pr} = \frac{\mu_\infty^* C_{p\infty}^*}{k_\infty^*}$,
and, hence,
\[
k_1 k_5 = \frac{\gamma}{\rm Pr_{local}}\ ,
\]
\noindent where 
\[
{\rm Pr_{local}} = \frac{\mu^* C_p^*}{k^*}\ .
\]

For a very wide range of temperatures, ${\rm Pr} \simeq 0.72$ (for air).
Since ${\rm Pr}^{-1} \simeq 1.4$, then the values
of $k_1$ and $k_5$ specified by the NAS benchmark are appropriate.

On the other hand, for most flows of interest, $\rm Re \sim O\left( 10^6
\right)$ or higher.  Hence, the relatively high values of $k_3 k_4$ does not
correspond well to the values in a typical ``real'' CFD code.  However, since
second-order central differences are used to discretize the governing
equations, a low value of $k_3 k_4$ would produce non-physical oscillatory
solutions.  Therefore, a relatively high value of $k_3 k_4$ is chosen to
circumvent such numerical difficulties in the solution of the synthetic
PDE's.
\newpage

\noindent{\bf (B)\ \ Additional terms in the NAS benchmark}
\vspace{10pt}

All such terms are associated with the constants $d_\xi^{\,(m)}$,
$d_\eta^{\,(m)}$ and $d_\zeta^{\,(m)}$, $m=1,5$.  An exact correspondence
with the original N-S equations would be possible if the above
constants are all zero; however, they are O(1) in the NAS benchmark.
Upon inspection of these terms, it is clear that they all represent an
added second-order dissipation although this has not been explicitly
stated in the text of the NAS-BM.  In practice, a blend of second and
fourth-order dissipation functions are usually incorporated into the
discretized equations to suppress high-frequency, non-physical
oscillations.  However, such a second-order dissipation function with
constants of O(1) would produce incorrect (excessively diffused)
solutions in a real CFD application.  Note, however, that the inclusion
of these terms does not alter the computational algorithm to solve
these equations nor does it alter the communication pattern.

\end{document}