%\documentstyle[11pt]{article} %\input{rep12} %\input{symbollga} %\begin{document} %\setcounter{page}{2} \newpage \begin{center} \Large{\bf II. RADIATION HYDRODYNAMICS} \end{center} \addcontentsline{toc}{section}{\protect\numberline{II.}{\bf ~Radiation Hydrodynamics}} \vspace*{5mm} {\flushleft\large{\bf A. The Equations}}\\ \vspace*{2mm} \addcontentsline{toc}{subsection}{\protect\numberline{II.A.}{~The Equations}} \begin{center} \it Reference \end{center} {\flushleft \hspace{10mm}\rm D. Mihalas and B. W. Mihalas, ``Foundations of Radiation\\ \hspace{15mm}Hydrodynamics'', Oxford University Press, New York, 1984}\\ \vspace{2mm} \noindent\rm The equations of Newtonian gray radiation hydrodynamics for a single fluid in the comoving-frame of reference taking into account terms of $O\Bigl({v\over c}\Bigr)$ are formulated according to Mihalas \& Mihalas. The two meaningful geometries for radiation transfer are slab symmetry, designated by the parameter $\mu = 0$, and spherical symmetry, $\mu = 2$. The hydrodynamical equations by themselves can also be solved in cylindrical symmetry, $\mu = 1$.\vspace*{2mm} \noindent\rm The mass within a shell of radius $r$ (a slab of distance $r$) is given by by the density: $$ m(r) = \int^r \rho(x) {d(x^{\mu+1})\over{\mu+1}} \eqno{(1)} $$ where $\rho$ is the density of the gas. The set of governing equations are, defining the comoving time derivative as $\Dt\equiv {\p\overp t} + u{\p\overp r}$: \begin{description} \item~{\sl continuity} $$ \Dt(\rho) + \rho {\p(r^\mu u)\over r^\mu\p r} = 0 \eqno{(2)} $$ \item~{\sl momentum} $$ \rho\Dt(u) + {\p(\Pg)\overp r} + {\p\Bigl(r^{3\mu/2}\PQ\Bigr)\over{r^{\mu/2}\p r}} + \rho \Bigl[\Bigl(1-{\mu\over2}\Bigr) g + 2\pi\mu {G m\over r^\mu}\Bigr] - {1\over c} \rho \kF \Fr = 0 \eqno{(3)} $$ \item~{\sl radiative momentum} $$ \rho\Dt\Bigl({\Fr\over{\rho c^2}}\Bigr) + \Bigl[{\p (\pr)\overp r} + {\mu\over2}{(3\pr - \er)\over r}\Bigr] + {\Fr\over c^2} {\p u\overp r} + {1\over c} \rho \kF \Fr = 0 \eqno{(4)} $$ \item~{\sl fluid plus radiative energy} $$ \rho\Dt\Bigl(\eg+{\er\over\rho}\Bigr) + {\p (r^\mu \Fr)\over{r^\mu\p r}} + (\Pg + \pr){\p(r^\mu u)\over r^\mu\p r} + \PQ \Bigl[{\p u \overp r}-{\mu\over2}{u\over r}\Bigr] $$ $$ + {\mu\over2}(\er - 3\pr){u\over r} = 0 \eqno{(5)} $$ \item~{\sl radiative energy} $$ \rho\Dt\Bigl({\er\over\rho}\Bigr) + {\p (r^\mu \Fr)\over{r^\mu\p r}} + \pr{\p(r^\mu u)\over r^\mu\p r} + {\mu\over2}(\er - 3\pr){u\over r} $$ $$ + c \rho \Bigl[\kE\er-\kP\ar T^4\Bigr] = 0 \eqno{(6)} $$ \end{description} \noindent\rm In equation (5) $\eg$, the specific internal energy of the gas, and the radiative energy $\er$, are combined to form a ``total internal'' energy. Correspondingly there are two different contributions to the work: that of the gas pressure~$\Pg$ and that of the radiative stress $\pr {\p(r^\mu u)\over r^\mu\p r} + {\mu\over2}(\er - 3\pr){u\over r}$. In addition we provide for the work done by an artificial stress $\PQ$. The radiant energy is transported via the radiative flux $\Fr$. We have neglected terms that combine radiation quantities with the gas acceleration.\vspace*{2mm} \noindent\rm The gas momentum equation (3) allows for gravitational forces in two geometries, with $G$ standing for the Newtonian gravitational constant, and $g$ for a constant acceleration. The radiation force, ${1\over c}\rho\kF\Fr$, the gradient of the gas pressure, and the artificial stress are exerted on the fluid elements, which move with velocity $u$.\vspace*{2mm} \noindent\rm The equation for the radiative flux (4) can be interpreted as the photon momentum equation which couples to the gas momentum via the last term.\vspace*{2mm} \noindent\rm The radiative energy equation (6) describes transport and work of the photons and constitutes, together with equation (4), the first two frequency and solid angle integrated moments of the transport equation. The last term in (6) can be viewed approximately as ${1\over {\tau_r}} \Bigl(\er - \ar T^4\Bigr)$, $\ar$ denoting the radiation constant, and then be interpreted as a heat exchange on a time scale $\tau_r \equiv {1\over {c\:\kappa\rho}}$, which is the lifetime of a photon traveling a mean free path ${1\over{\kappa\rho}}$ between two scattering or absorption events. If $\er > \ar T^4$ then radiative energy is converted into heat, and {\it vice versa}. By assuming such a form for the radiative coupling term we restrict ourselves to a radiating fluid in LTE.\vspace*{2mm} \noindent\rm Equations (4) \& (6) are referred as ``full transport'' which is a correct description in all regimes for the one-dimensional case. In particular, the ``diffusion regime'' is recovered when (4) becomes equivalent to $$ \Fr = - {c \over 3\rho\chi}~{\p\er\overp r}. \eqno{(7)} $$ Substituting formula (7) into equation (6) yields the non-equilibrium diffusion because $\ar T_r^4 \equiv \er$ defines a radiative temperature field which is allowed to differ from the fluid temperature $T$. Imposing the additional condition that $\kE\er = \kP\ar T^4$, or $T_r = T$ leads back to equilibrium diffusion.\vspace*{2mm} \noindent\rm In spherical symmetry it is customary to define the luminosity as the surface brightness ${\cal L} = 4\pi r^2 \Fr$ of a shell with radius $r$. \vspace*{5mm} {\flushleft \large{\bf B. Adaptive Mesh Transformation}}\\ \vspace*{2mm} \addcontentsline{toc}{subsection}{\protect\numberline{II.B.}{~Adaptive Mesh Transformation}} \begin{center} \it Reference \end{center} {\flushleft \hspace{10mm}\rm Adaptive--Mesh Radiation Hydrodynamics. I. The Radiation\\ \hspace{15mm}Transport Equation in a Completely Adaptive Coordinate System,\\ \hspace{15mm}K.--H. Winkler, M.L. Norman, and D. Mihalas, \it J.Q.S.R.T.,\\ \hspace{15mm}\bf 31, \rm 473, 1984\\ \hspace{10mm}\rm W. M. Tscharnuter and K.--H. A. Winkler, \it Comput. Phys. \\ \hspace{15mm}Commun., \bf 18, \rm 171, 1979\\} \vspace*{2mm} \noindent\rm The mass definition (1) can be immediately rewritten in differential form as: $$ dm = \rho~d\Bigl[{r^{\mu+1}\over{\mu+1}}\Bigr] \eqno{(8)} $$\vspace*{2mm} \noindent\rm The partial differential equations (2)--(6) should be represented in finite difference form on an arbitrary moving grid distribution. This distribution can be obtained by uniquely subdividing the whole computational domain into cells around each grid point $k$ with $k = 1, \dots N$. If (2)--(6) are integrated over those cells each encompassing a volume (length) $dV_k = {r^{\mu+1}\over{\mu+1}}\Big|_k$, they reduce to a set of $N$ ordinary differential equations (in time). The proper transformation is achieved by the adaptive--mesh Reynolds transport theorem as described in Winkler~\etal Hereby the cells are assumed to move with a relative velocity $\ur|_k$ with respect to the fluid flow.\vspace*{2mm} \noindent\rm For each cell $(k = 1, \dots N)$ the following equations hold: \begin{description} \item~{\sl continuity} $$ {d\over d t} \int_V \rho~dV + \oint_{\p V}~\rho\ur~dS = 0 \eqno{(9)} $$ \item~{\sl fluid momentum} $$ {d\over d t} \int_V u \rho~dV - \oint_{\p V}~u \rho\ur~dS + \int_V r^\mu d(\Pg) + \int_V r^{-\mu/2} d[r^{3\mu/2}~\PQ] $$ $$ + \int_V \Bigl[\Bigl(1-{\mu\over2}\Bigr) g + 2\pi\mu {G m\over r^\mu} - \kF {\Fr\over c}\Bigr]~\rho dV = 0 \eqno{(10)} $$ \item~{\sl radiative momentum} $$ {d\over d t} \int_V {\Fr\over c^2}~dV - \oint_{\p V}~{\Fr\over c^2}\ur~dS + \int_V r^\mu d(\pr) + \int_V~\Bigl[{\mu\over2}{(3\pr-\er)\over r}\Bigr]~dV $$ $$ + \int_V {\Fr\over c^2} r^\mu d(u) + \int_V \Bigl[\kF {\Fr\over c}\Bigr]~\rho dV = 0 \eqno{(11)} $$ \item~{\sl fluid plus radiative energy} $$ {d\over d t} \int_V \Bigl(\eg+{\er\over\rho}\Bigr) \rho~dV - \oint_{\p V}~\Bigl(\eg+{\er\over\rho}\Bigr)\rho\ur~dS + \int_V d[r^\mu~\Fr] $$ $$ + \int_V (\Pg + \pr)~d[r^\mu~u] + \int_V \PQ\Bigl[{\p u\overp r}-{\mu\over2}{u\over r}\Bigr]~dV + \int_V \Bigl[{\mu\over2}(\er-3\pr){u\over r}\Bigr]~dV = 0 \eqno{(12)} $$ \item~{\sl radiative energy} $$ {d\over d t} \int_V \er ~dV - \oint_{\p V}~\er \ur~dS + \int_V d[r^\mu~\Fr] $$ $$ + \int_V \pr~d[r^\mu~u] + \int_V \Bigl[{\mu\over2}(\er-3\pr){u\over r}\Bigr]~dV + \int_V~c~\Bigl[\kE\er-\kP\ar T^4\Bigr]~\rho dV = 0 \eqno{(13)} $$ \end{description} \noindent\rm Here $\p V$ denotes the surface of the volume element $V$. In the transformation onto the adaptive grid Neumann-Richtmyer's idea of introducing an artificial viscosity to treat shocks in the flow has been implemented according to Tscharnuter \& Winkler. A tensor formulation is used which introduces the artificial stress: $$ \PQ = - {4\over3} \rho \mq \Bigl[{\p u\overp r} - {\mu\over2}{u\over r}\Bigr] \eqno{(14)} $$ in conjunction with a linear and quadratic artificial viscosity, $\ellQ$ giving the shock width: $$ \mq = C_{\scriptscriptstyle 1} \ellQ \sqrt{\gamma\Pg/\rho} - C_{\scriptscriptstyle 2} \ellQ^2 \min\Bigl[~0, {\p(r^\mu u)\over{r^\mu\p r}}~\Bigr] \eqno{(15)} $$ \vspace*{5mm} {\flushleft \large{\bf C. Symbolic Notation}}\\ \vspace*{2mm} \addcontentsline{toc}{subsection}{\protect\numberline{II.C.}{~Symbolic Notation}} \begin{center} \it Reference \end{center} {\flushleft \hspace{10mm}\rm WH80s: Numerical Radiation Hydrodynamics, K.--H. Winkler and \\ \hspace{15mm}M.L. Norman,in \it Astrophysical Radiation Hydrodynamics,\\ \hspace{15mm}\rm (Dordrecht: Reidel), 71, 1986}\\ \vspace*{2mm} \noindent\rm The volume (length) integrals in the equations above (9)--(13) can be approximated with finite differences in the simplest way if one assumes that the integrands remain constant. In general this is valid as long as the linear dimensions of the integral volumes are small compared to any physical length scale in the problem. In the case of smooth variations of the physical variables this can be achieved by decomposing the physical domain into a sufficiently large number of numerical zones. In the case of ({\em in praxi}) discontinuities the demand is that the physical variables still are not allowed to vary significantly between neighboring zones. This can be satisfied only if the zones become very small; or in other words if the numerical grid {\em resolves} such nonlinear features. Both situations are taken care of by the implicit nature of the adaptive mesh.\vspace*{2mm} \noindent\rm In a like spirit the spatial differential operators can interpreted in the simplest form by spatial differences denoted by $\Delta$ in symbolic form. Equation (8) valid for all zones $k = 1, \dots, N$ on the numerical domain then can be straightforwardly expressed as: $$ \Delta m = \rho \Delta V \eqno{(16)} $$\vspace*{2mm} \noindent\rm The temporal derivatives remain intact, the surface integrals are resolved as differences of the surface terms, and the spatial derivatives approximated as spatial differences the symbolic form of the equations (9)--(15) are again valid for all zones $k = 1, \dots, N$: \begin{description} \item~{\sl continuity} $$ \frac{d(\rho\Delta V)}{dt} + \Delta(\rho r^{\mu} \ur) = \Delta\left(r^{\mu}\sr \frac{\Delta\rho}{\Delta r}\right) \eqno{(17)} $$ \item~{\sl fluid momentum} $$ \frac{d}{dt}[u\Big<\rho\Delta V\Big>] - \Delta\left(\Big<\frac{dm}{dt}\Big> u\right) + r^{\mu}\Delta\Pg + \left(\frac{2-\mu}{2}g + \frac{2\pi\mu G m}{r^{\mu}}\right) \Big<\rho\Delta V\Big> $$ $$ = \phi_{Q}\Delta V + \frac{1}{c}<\!\!\kF\!\!>\Fr\Big<\rho\Delta V\Big> \eqno{(18)} $$ \item~{\sl radiative momentum} $$ \frac{d}{dt} \left[ \frac{\Fr}{c^2}\Big<\Delta V\Big> \right] - \Delta\left[\Big<\frac{dm}{dt}\Big>\left(\frac{\Fr}{c^2<\!\!\rho\!\!>}\right)\right] + r^{\mu}\left(\Delta\pr + \frac{\Fr\Delta u}{c^2}\right) + \frac{(3\pr - \er)}{r}\Big<\Delta V\Big> $$ $$ = - \frac{1}{c}<\!\!\kF\!\!>\Fr\Big<\rho\Delta V\Big> \eqno{(19)} $$ \item~{\sl fluid plus radiative energy} $$ \frac{d}{dt}\left[\left(\eg +\frac{\er}{\rho}\right)\rho\Delta V\right] - \Delta\left[\frac{dm}{dt}\left(\eg +\frac{\er}{\rho}\right) - r^{\mu}\Fr\right] $$ $$ + (\Pg+\pr)\Delta(r^{\mu}u) + (\er-3\pr)\frac{u}{r}\Delta V = \Delta\left(r^{\mu}\se\frac{\Delta\eg}{\Delta r}\right) + \epsilon_{Q}\Delta V \eqno{(20)} $$ \item~{\sl radiative energy} $$ \frac{d}{dt} \left( \er \Delta V \right) - \Delta\left[\frac{dm}{dt}\left(\frac{\er}{\rho}\right) - r^{\mu}\Fr \right] + \pr\Delta(r^{\mu}u) + (\er - 3\pr)\frac{u}{r}\Delta V $$ $$ = c \Bigl[\kP\ar T^4 - \kE\er\Bigr]\rho\Delta V \eqno{(21)} $$ \end{description} \noindent\rm The bracketed quantities are averaged in this staggered mesh representation, \cf~the discussion in section III A. This set of equations corresponds to the equations $M1$, $C1$, $GM1$, $RM1$, $FE1$, $RE1$ in the TITAN Code Reference Manual.\vspace*{2mm} \noindent\rm The original advection operator $\Delta(\rho r^{\mu} \ur)$ has been retained only in the continuity equation (17). In all other equations, for reasons discussed by Winkler and Norman, the combination $\rho r^{\mu} \ur$ is replaced by $- \frac{dm}{dt}$. This is possible because the motion of the grid relative to the fluid results in an apparent mass transfer rate that can be computed according to: $$ \frac{dm}{dt} = - \rho r^{\mu} \ur \eqno{(22)} $$ Since the correct amount of mass per zone is already determined by (16) the advection operator in all other ordinary differentail equations (17)--(21) can be expressed through $-\Delta\left(\frac{dm}{dt}\right)$.\vspace*{2mm} \noindent\rm The artificial stress enters into the momentum equation (18) in the form of a force: $$ \phi_Q \Delta V \equiv r^{-(\mu/2)} \Delta\left[\frac{4}{3} \rho \mq r^{(3\mu/2)} \left(\frac{\Delta u}{\Delta r} - \frac{\mu}{2}\Bigl<\frac{u}{r}\Bigr>\right)\right] \eqno{(23)} $$ and into the fluid plus radiative energy equation (20) in the form of a source: $$ \epsilon_Q \equiv \frac{4}{3}\mq\rho\left(\frac{\Delta u}{\Delta r} - \frac{\mu}{2}\Bigl<\frac{u}{r}\Bigr>\right)^{2} \eqno{(24)} $$ The artificial viscosity coefficient is evaluated according to: $$ \mq \equiv C_{\scriptscriptstyle 1} \ellQ a_s - \min [ C_{\scriptscriptstyle 2} \ellQ^2 \Delta(r^{\mu}u) / \Delta V, 0]\eqno{(25)} $$ The shock width can be assumed as a constant $\ellQ = \ell_0$ or proportional to the radius $\ellQ = \ell_1 <\!\!r\!\!>$. Additional diffusion operators appear in the continuity and fluid plus radiative energy equations. Their purpose is to treat numerical artifacts that occur in contact continuities and when shock reflect from walls (wall heating). The constants $\sr$ and $\se$ control the length scale over which density or internal energy is allowed to diffuse. \vspace*{5mm} {\flushleft \large{\bf D. Boundary Conditions}}\\ \vspace*{2mm} \addcontentsline{toc}{subsection}{\protect\numberline{II.D.}{~Boundary Conditions}} \noindent\rm TITAN is equipped with a broad variety of boundary conditions, allowing for everything but very special-purpose boundaries. They naturally fall into two categories which describe the behavior of the radiation field and of the fluid ath the boundaries. The inner and outer (left and right) boundary conditions obey similar laws; hence is suffices to give just the generic formulation. More details are given in the TITAN Code Reference Manual.\\ \noindent (1) radiative boundary:\vspace{2mm} \noindent In general the radiative boundary can be posed by requiring the radiative energy to remain constant: $\Delta\er = 0$. We provide for the optically transmitting and reflecting cases. The first is defined through a flux which is computed from the free streaming formula. In addition the flux is allowed to be modified by a constant intensity field $I^\pm$. The second case is given through a zero flux condition. Furthermore we have the possibility of imposing a time dependent net flux. \begin{description} \item~{\sl transmitting} $$ \Fr = \pm c \gE \er \mp (2 \gE + 1) \pi I^\pm. \eqno{(26)} $$ \item~{\sl reflecting~flux} $$ \Fr = 0. \eqno{(27)} $$ \item~{\sl imposed~flux} $$ \Fr = \Phi(t). \eqno{(28)} $$ \end{description} \noindent In equation (26) $\gE$ is a surface Eddington factor connecting $\er$ to $\Fr$.\\ \noindent (2) Eulerian boundary:\vspace{2mm} \noindent The Eulerian boundary is given by the fixed grid condition $\dot r = 0$. We allow for three different cases. The transmitting condition demands that temperature, density, and velocity all remain constant across the boundary. The zero flux condition necessitates that the temperature and density stay constant while the velocity vanishes. For the nonzero flux condition density and temperature can be functions of time. \begin{description} \item~{\sl transmitting} $$ \Delta T = 0~;~\Delta\rho = 0~;~\Delta u = 0. \eqno{(29)} $$ \item~{\sl zero~flux} $$ \Delta T = 0~;~\Delta\rho = 0~;~u = 0. \eqno{(30)} $$ \item~{\sl nonzero~flux} $$ T = \Theta(t)~;~\rho = P(t)~;~u = 0. \eqno{(31)} $$ \end{description} \noindent (3) Lagrangean boundary:\vspace{2mm} \noindent The Lagrangean boundary is given by the fixed mass condition $\dot m = 0$. We allow for two different forcings. In the piston driven case a time dependent velocity is assumed to move the boundary. Otherwise we allow for a forcing through an external pressure. \begin{description} \item~{\sl external velocity} $$ u = U_{\scriptscriptstyle ext}(t). \eqno{(32)} $$ \item~{\sl external pressure} $$ \Pg = P_{\scriptscriptstyle ext}(t). \eqno{(33)} $$ \end{description} \noindent\rm In addition to the boundaries we need to provide for phantom zones so that the difference stencil stays well defined. These phantom zones have to be chosen to reflect the physical properties of the boundary conditions. \vspace*{5mm} {\flushleft \large{\bf E. Constitutive Relations}}\\ \vspace*{2mm} \addcontentsline{toc}{subsection}{\protect\numberline{II.E.}{~Constitutive Relations}} \begin{center} \it References \end{center} {\flushleft \hspace{10mm}\rm D. Mihalas, W. D\"appen, D. Hummer, and B. W. Mihalas,\\ \hspace{15mm}\it Ap. J., \bf 331, \rm 794, 1988; \it Ap. J., \bf 331, \rm 815, 1988;\\ \hspace{15mm}\it Ap. J., \bf 332, \rm 261, 1988; \it Ap. J., \bf 350, \rm 350, 1990}\\ \vspace{2mm} \noindent\rm (1) fluid equation of state:\vspace{2mm} \noindent\rm Arbitrary equations of state as functions of the density and gas temperature are permitted: \begin{description} \item~{\sl caloric} $$ \eg = \eg(\rho,T), \eqno{(34)} $$ \item~{\sl thermal} $$ \Pg = \Pg(\rho,T). \eqno{(35)} $$ \end{description} \noindent\rm They can be given in form of tables or as analytical formulae. In the first case monotonic bicubic hermite polynomials are employed to perform the interpolations. Important thermodynamical properties for the fluid are derived through the Gibbs form $d\eg(\rho,T) = T ds(\rho,T) + \Pg(\rho,T) d({1\over\rho})$, where $s$ denotes the specific entropy. Among those of particular interest is the adiabatic index $\gamma = {\partial \ln\Pg \overp {\ln\rho}}\Big|_s$.\vspace{2mm} \noindent\rm For the purpose of the test problems a perfect gas is assumed so that: $\Pg = {{\it R}\over\overline{\mu}}\rho T$, with ${\it R}$ standing for the universal gas constant and $\overline{\mu}$ for the mean molecular weight. Prescribing a constant adiabatic index $\gamma$ then relates the gas pressure to the gas energy: $\eg = {1\over (\gamma-1)}{\Pg\over\rho}$. \\ \noindent\rm (2) radiative equation of state:\vspace{2mm} \noindent The two-equation formalism for the radiative transport, equations (6) \& (4), is closed via the variable Eddington factor~$\fe$. The radiative pressure~$\pr$ in this description is simply proportional to $\er$: $$ \pr = \fe \er. \eqno{(36)} $$ In order to evaluate $\fe$ we integrate the static transport equation for the radiative intensity as a function of angle $\phi$, $I(x=\cos\phi)$, using a ray tracing technique. The Eddington factor is then computed as the ratio of the second to the zeroth moment: $$ \fe = {{\int_{-1}^{+1} I(x) x^2 dx}\over {\int_{-1}^{+1} I(x) dx}}. \eqno{(37)} $$ At the boundaries of the radiation field a corresponding factor, called $\gE$, is evaluated. $\fe$ describes the geometric properties of the radiation field. For an isotropic intensity distribution, $I(\phi) = I_o$, we obtain $\fe = {1\over3}$. This corresponds to the diffusion limit of the two-equation formalism. If the intensity distribution peaks along the line of sight, \ex~$I(\phi) \propto \cos^2\phi$, we get $\fe = {3\over5} > {1\over3}$. The stronger the intensity field is pointed along the line of sight the closer $\fe \rightarrow 1$. This corresponds to the free streaming limit of the full transport formalism. On the other hand, if the intensity distribution peaks sideways, \ex~ $I(\phi) \propto \sin^2\phi$, we find $\fe = {3\over15} < {1\over3}$. This situation is often found in radiating shock fronts. \\ \noindent\rm (3) opacities:\vspace{2mm} \noindent In general the frequency dependent opacity $\chi(\nu)$ is given by contributions from absorption, $\chi^a(\nu)$, and scattering, $\chi^s(\nu)$. The gray formulation of the radiation hydrodynamics equations introduces three different opacity means: the flux mean, the absorption (or energy) mean, and the Planck mean: \vspace*{2mm} \begin{description} \item {\sl flux~mean} $$ \kF = \int_{0}^{\infty} [\chi^a(\nu)+\chi^s(\nu)]~{\Fr(\nu)\over\Fr}~d\nu. \eqno{(38)} $$ \item {\sl absorption~mean} $$ \kE = \int_{0}^{\infty} \chi^a(\nu)~{\er(\nu)\over\er}~d\nu. \eqno{(39)} $$ \item {\sl Planck~mean} $$ \kP = \int_{0}^{\infty} \chi^a(\nu)~{B(\nu,T)\over \sigma T^4/\pi}~d\nu. \eqno{(40)} $$ \end{description} The flux and the energy mean both depend on the spectral distribution of the radiation field, whereas the Planck mean is given by the local gas temperature and density via the Planck distribution $B$. The evaluation of the first two opacity means therefore can be carried out correctly only if one also computes the frequency dependent, \ie~non-gray, radiation field. This is often done with a multigroup approach. For the moment TITAN is restricted to the gray transport. \vspace{2mm} \noindent\rm Defining the Rosseland mean opacity $$ \kR^{-1} = \int_{0}^{\infty} [\chi^a(\nu)+\chi^s(\nu)]^{-1}~ {{\p B(\nu,T)/\p T}\over 4\sigma T^3/\pi}~d\nu, \eqno{(41)} $$ one can show that in diffusion approximation it equals to the flux mean, \ie~$\kF = \kR$. For the purpose of our test problems we assume a constant $\kR$. In order to simulate the effects of a scattering or absorption dominated fluid we set both the energy mean and the Planck mean equal to a fraction $0 \le \fr \le 1$ of the Rosseland mean, \ie~$\kE = \kP = \fr \kR$.\\ \noindent\rm Both the fluid equation of state and the opacities are available in table lookup form. The source is described by Mihalas~\etal in a series of articles. In this case monotonized bicubic hermite interpolation is used to extract the desired quantity $x(T,\rho)$ along with its partial derivatives ${\p\ln x\overp \ln T}$ and ${\p\ln x\overp \ln\rho}$. %\end{document} .