\documentclass[11pt,english]{article}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
% \usepackage[cp850]{inputenc}
\usepackage{amsmath} % need for subequations
\usepackage{graphicx} % need for figures
% \usepackage{verbatim} % useful for program listings
% \usepackage{color} % use if color is used in text
% \usepackage{subfigure} % use for side-by-side figures
\usepackage{hyperref} % use for hypertext links, including those to external documents and URLs
\usepackage{rotating}
\usepackage{graphics}
\usepackage{amssymb}
\usepackage{listings} \lstset{language=Fortran}
\usepackage{fullpage}
\usepackage{wasysym}
% \usepackage{pstricks}
\usepackage{color}
\usepackage{float}
\usepackage{helvet}
\makeatletter
\newcommand{\rmnum}[1]{\romannumeral #1}
\newcommand{\Rmnum}[1]{\expandafter\@slowromancap\romannumeral #1@}
\makeatother
% **********************************************
% Layout:
% http://en.wikibooks.org/wiki/LaTeX/Page_Layout
% **********************************************
\setlength{\parindent}{0pt}
\setlength{\parskip}{8pt}
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}
\begin{document}
\thispagestyle{empty}
\begin{center}
% \includegraphics[width=\textwidth]{SN.Cover.eps}
\includegraphics[width=16.2cm]{SN.Cover.eps}
\end{center}
\newpage
\thispagestyle{empty}
.
\cleardoublepage
\newpage
\title{\bf Numerical Solutions of Partial Differential Equations Using a Hydrodynamic Supernova Model}
\author{\bf Michael Kenn}
\date\today
\maketitle
\begin{abstract}
This thesis deals with how to solve nonlinear partial differential equations numerically.
For a supernova model, I examine the compressible hydrodynamic Euler equations with non-continuous initial conditions.
I compare highly sophisticated differential schemes and combine them with smoothing techniques to repress oscillations.
The computational results are compared with theoretical models where possible.
After having established a reliable framework, the effects of different initial conditions, heat conduction, artificial viscosity and the like are studied.
% , to estimate the reliability of the results.
% Finally I use some of the best approaches to simulate astrophysical questions such as the collision of supernova remnants in star generation regions.
\end{abstract}
\section*{Preface}\nonumber
In recent years techniques using numerical schemes for solving nonlinear hyperbolic PDEs have improved vastly.
There are many different approaches to repressing oscillations, which are a direct consequence of using higher-order schemes.
Problems including non-continuous initial conditions, inclusion of source terms due to cylindrical or spherical symmetry or physical effects such as heat conduction or viscosity, can be attacked extensively with a pool of appropriate methods.
However, not all methods are suitable for certain models, for example a supernova simulation, and it is not easy to know which method may be the best to start with.
For the unacquainted student, dealing with the many possible techniques can be irritating at a first glance.
Even though there is a large amount of literature that deals with all kinds of special cases, it is hard to decide which approach will be the most
% prosperous
successful for a particular problem.
I found myself in this position when I first started investigating various differential schemes to solve compressible hydrodynamic Euler equations (a set of highly nonlinear PDEs) to simulate the expansion of supernova remnants.
While studying the literature, I found some prescribed solutions \cite{CIO001}.
Had it only been only about producing a hydrodynamic code to run some simulations,
% one could easily take the present suggestion as a black box.
one could have taken the easier way, i.e.\ taking the present suggestion as a black box.
However to acquire a deeper understanding of the matter I knew I had to take the hard way.
I decided to understand why certain approaches are superior than others and why some methods fail.
The underlying objective of this thesis therefore is to develop an in-depth insight into numerical simulations of flows in astrophysics rather then to use a pre-defined solution.
The goal of this thesis is to give a comprehensive summary of all issues required for an understanding of a simple supernova simulation.
Chapter \ref{SEC001} gives a short overview of my strategy.
Chapter \ref{SEC002} is a crash course in astrophysics, covering the Euler equations, the expansion of supernova remnants (SNR), shock waves and the Rankine-Hugoniot conditions, the composition of the interstellar medium (ISM) and so on.
Chapter \ref{SEC003} gives an introduction to the theory of numerical solutions to hyperbolic PDEs.
I will introduce a large number of differential schemes and inspect them on a sample-model, where exact solutions are available in chapter \ref{SEC004}.
% Chapter \ref{SEC004} is a test case for the numerical schemes introduced in \ref{SEC003}.
Chapter \ref{SEC005} finally deals with the supernova model.
Features such as heat conduction, viscosity and turbulences are investigated and compared to predictions based on astrophysical considerations.
Chapter \ref{SEC006} gives prospects for possible further investigations.
% With the certainty of having chosen a suitable approach in chapter \ref{SEC005}, chapter \ref{SEC006} will simulate the agglomeration caused by the collision of SNRs, which is required for further star generation.
Appendix \ref{SEC007} gives a brief program documentation.
% At this point I would like to express my gratitude to two persons who made this work possible.
At this point I would like to express my gratitude towards two persons for their extensive and valuable support during the entire time I worked on this thesis.
I am extremely grateful to Prof. Neumann of the Department of Computational Physics at the University of Vienna, who took the time for weekly meetings with me to work through complex issues and always delivered the right ideas whenever I did not know the way forward.
The other person is my wife Ritu, who not only made it possible for me to return to academics and pursue my passion, but also spent a considerable amount of time to proof-read my thesis.
% but also spent a considerable amount of her time supporting me to complete my thesis.
% and was also instrumental in the completion of this thesis.
\vspace{20mm}
\begin{flushright}
Michael Kenn\\
Vienna, May 2014
\end{flushright}
\newpage
\section{Introduction}\label{SEC001}
\renewcommand{\abstractname}{}
\begin{abstract}
{\it
This chapter explains the reason for the chosen way to develop a suitable framework for astrophysical simulations.
The single steps are introduced individually and dependencies are identified.
There is a strong focus on variety and reliability of the applied methods rather than using a ready-to-use solution.
}
\end{abstract}
\subsection{Workflow}
We should be aware that the final goal is to develop a framework to solve a set of (a) highly nonlinear PDEs with (b) possibly non-continuous initial conditions under (c) difficult boundary conditions using (d) certain geometric constraints (e.g.\ spherical symmetry).
It is therefore reasonable to break down the whole issue into small steps.
\subsubsection{Physical issues}
First of all it is necessary to identify the underlying physics.
In the case of a supernova explosion, as we will see in chapter \ref{SEC002}, it is a good approach to distinguish between the following states:
\paragraph{Initial state:}
We assume a self-similar expansion in a short initial phase (e.g.\ 200 years), for which we need to chose a fixed initial density profile.
An initial value problem with constant initial conditions except one discontinuity is called a {\it Riemann problem} (RP).
In the case of the Euler equations in one linear dimension the Riemann problem can be solved analytically (e.g.\ \cite{TOR001}).
We will take advantage of this and check our numerical schemes with these well known solutions.
Only then we will expand to spherical symmetry and choose particular, non-constant initial profiles for density, pressure and energy.
\paragraph{Adiabatic expansion:}
The first stage of the simulation is purely pressure driven and adiabatic.
No heat conduction or energy loss has to be considered.
The lack of a cooling function simplifies the set of equations considerably.
Theoretical predictions can be used for checking the quality of the model.
At this point we will introduce the concept of artificial viscosity and adjust the model parameters to match the theoretical predictions.
\paragraph{Heat conduction:}
As we have a simple but accurate model, we will now introduce perturbation terms.
This will primarily be a cooling function modeling the effect of heat conduction and viscosity.
The Euler equations for momentum conservation then turn into the Navier--Stokes equations (\cite{TOR001}, p17), which are, as is well-known, hard to tackle.
\subsubsection{Numerical issues}
There are as many different optimal solutions as there are problems.
Even when it becomes clear which type of differential scheme might be the best fit, parameter tuning can still be very cumbersome in a high-dimensional parameter space.
To find a suitable scheme my approach is simple elimination.
I start with approximately 15 differential schemes, apply them to the Euler equations and challenge them by continuously adding new requirements.
There are a few major hurdles to master:
\paragraph{Riemann problem without boundary conditions:}
The differential scheme must, of course, solve the one-dimensional, linear Riemann problem with sufficient accuracy.
First order schemes such as the Godunov scheme are sorted out at this point.
\paragraph{Spherical symmetry:}
Spherical (and cylindrical) symmetry can be implemented by adding a source term (this will be discussed in chapter \ref{SEC003}).
Some differential schemes are very sensitive towards the origin and therefore inapplicable for supernova simulations.
In any case, much care must be taken while handling the boundary at the origin.
\paragraph{Supernova features:}
A major problem in the implementation of the supernova model is a large range of magnitudes for the parameters.
The radial velocity, especially at the beginning, is very high and there is a tremendous increase in pressure at the expanding shell.
Numerical oscillations can easily lead to negative pressure or negative internal energy, which is not only non-physical, but will in most cases also cause the divergence of the scheme.
Second order schemes not taking advantage of total variation diminishing (TVD) techniques will not pass this point.
\paragraph{Extensions of the model:}
Finally heat conduction, viscosity (and possibly turbulence) has to be considered to provide a good model.
Not all of the remaining schemes can deal with these requirements in a sufficient way and here a final choice is made.
\subsubsection{Scientific issues}
The main objective of this work is to choose an appropriate framework to be used for further investigations.
The material at hand should be a stepping stone to give up spherical symmetry and switch to a comprehensive three-dimensional model.
With this extended model physical processes which could not be covered here can be studied.
% Once we settle for an appropriate model, we use it for further investigations.
% At this point we give up spherical symmetry and switch to a comprehensive three-dimensional model.
% We verify that the three-dimensional model actually is an extension of the original model, i.e.\ if the results can be reproduced from the spherically symmetrical case.
% With this model, the collision of supernova waves can be studied.
% These are some of the reasons for matter densities in star-forming regions responsible for new star generations.
\newpage
\section{Astrophysical background}\label{SEC002}
\renewcommand{\abstractname}{}
\begin{abstract}
{\it
This chapter gives an overview of the astrophysics required to understand the expansion of supernova remnants (SNR).
I will begin with the Euler equations for compressible fluids.
After that I will discuss properties of shock waves, including the Rankine-Hugoniot conditions, which describe the contact discontinuity caused by an infinitesimally thin shock wave.
There will be an introduction to the interstellar medium (ISM), covering the heating and cooling processes, which are essential for the energy loss of the SNR.
A major part will be the different phases of the expansion of the SNR.
These are very important as here we obtain theoretical information which should match with the numerical model.
}
\end{abstract}
\subsection{Euler equations}\label{EULER}
Depending on the author's area of expertise, the Euler equations may be written in a number of ways.
A classical approach is that of Landau \& Lifschitz \cite{LAN001}, where the Euler equations describe the behaviour of ideal fluids.
In astrophysics, the Euler equations characterize the motion and conservation of matter as outlined by Schneider \cite{SCN001}.
A problem-specific transformation of the Euler equations is used in the paper of Sedov \cite{SED001}, which is advantageous for the derivation of an exact blast solution.
Last but not least, a formulation useful for numerical calculations is given by Toro \cite{TOR001}.
This chapter attempts to compile all approaches above into a comprehensive description.
\subsubsection{Formulation of the Euler equations}\label{EULERSUB}
The Euler equations are a system of hyperbolic partial differential equations to describe compressible, inviscid flows.
One popular way to formulate them is in {\it conservative} form
\begin{eqnarray}
\partial_t\rho+\vec\nabla(\rho\mathbf v)&=&0\label{EQU001}\\
\partial_t(\rho\mathbf v)+\vec\nabla((\rho\mathbf v)\mathbf v+p\mathbf I)&=&\mathbf{0}\label{EQU002}\\
\partial_t\epsilon+\vec\nabla(\epsilon\mathbf v+p\mathbf v)&=&0\label{EQU003}
\end{eqnarray}
where
\begin{itemize}
\item $\rho$ is the mass density of the fluid
\item $\mathbf v=(v_1,v_2,v_3)$ is the velocity vector of the fluid
\item $\epsilon$ is the total specific energy density of the fluid
\item $p$ is the pressure
\end{itemize}
The term {\it conservative} implies that the equations (\ref{EQU001})-(\ref{EQU003}) written in integral form describe conservation of mass, momentum and energy.
The notation of energy is not unique in the literature and here I will use the following conventions to avoid ambiguities:
\begin{center}
\begin{tabular}{r|l}
total energy density &
$\epsilon=\rho(\frac{1}{2}\mathbf v^2+U)=\rho E$ \\
internal energy density &
$u=\rho U$ \\
total specific energy &
$E=\frac{1}{2}\mathbf v^2+U=T+U$ \\
specific internal energy &
$U=\frac{\epsilon}{\rho}-\frac{1}{2}\mathbf v^2$ \\
specific kinetic energy &
$T=\frac{1}{2}\mathbf v^2$ \\
\end{tabular}
\end{center}
It is important to understand that equations (\ref{EQU001})-(\ref{EQU003}) form a set of $1+n+1=5$ equations in the six variables $\rho, v_1, v_2, v_3, \epsilon$ and $p$, where $n=3$ is the number of spatial dimensions.
An equation of state (EOS) $p=p(\rho,\epsilon)$ which correlates the pressure $p$ to the specific volume $v=\frac{1}{\rho}$ und the total energy density $\epsilon$ is required in addition.
The ISM is a very thin medium and we can therefore assume the EOS for an ideal gas
% \footnote{The case of a co-volume gas can be treated similarly, the vacuum case needs some additional care \cite{TOR001}.}
% co-volume gas: Toro p36 (1.50)
\begin{eqnarray}
p &=& (\gamma-1)\rho U=(\gamma-1)(\epsilon-\frac{1}{2}\rho\mathbf v^2)\label{EQU004}\\
\epsilon &=&\frac{p}{\gamma-1}+\frac{1}{2}\rho\mathbf v^2
\end{eqnarray}
with $\gamma$ being the adiabatic index (heat capacity ratio).
The Euler equations can also be formulated in vector form as
\begin{eqnarray}
\partial_t\mathbf u+\vec\nabla f(\mathbf u)=0,\qquad
\mathbf u=\begin{pmatrix}\rho\\\rho\mathbf v\\\epsilon\end{pmatrix},\qquad
f(\mathbf u)=\begin{pmatrix}\rho\mathbf v\\
(\rho\mathbf v)\mathbf v+p\mathbf I\\
(\epsilon+p)\mathbf v
\end{pmatrix}\label{EQU005}
\end{eqnarray}
Equation (\ref{EQU005}) represents the conservative form of the Euler equations and we will refer to the set $\mathbf u=(\rho,\rho\mathbf v,\epsilon)$ as the conservable variables with respect to mass, momentum and energy.
Another alternative formulation of the Euler equations -- under the assumption of an ideal gas -- is given by
\begin{eqnarray}
\partial_t\mathbf w+F(\mathbf w)\vec\nabla\mathbf w=0,\qquad
\mathbf w=\begin{pmatrix}\rho\\\mathbf v\\p\end{pmatrix},\qquad
F(\mathbf w)=\begin{pmatrix}
\mathbf v & \rho & 0 \\
0 & \mathbf v & \frac{1}{\rho} \\
0 & a^2\rho & \mathbf v
\end{pmatrix}\label{EQU006}
\end{eqnarray}
where $a=\sqrt{\frac{\gamma p}{\rho}}$ is the adiabatic sound speed.
The second component of (\ref{EQU006}) should be read as
\begin{eqnarray}
\partial_t\mathbf v+(\mathbf v\cdot\vec\nabla)\mathbf v&=&-\frac{\nabla p}{\rho}
\end{eqnarray}
and the third component of (\ref{EQU006}) together with (\ref{EQU004}) can be also written as
\begin{eqnarray}
\partial_t(\rho U)+\vec\nabla(\rho U\mathbf v)&=&-p\vec\nabla\mathbf v
\end{eqnarray}
Equation (\ref{EQU006}) is called the {\it non-conservative} form of the Euler equations, and the set $\mathbf w=(\rho,\mathbf v,p)$ is called the {\it primitive} variables.
Another formulation of the third equation under the assumption of an ideal gas
% ($p/\rho=(\gamma-1)U$)
which is common in the literature \cite{SED001}\cite{SED002}\cite{KAM001} is
\begin{eqnarray}
\partial_t(\frac{p}{\rho^\gamma})+\mathbf v\cdot\vec\nabla(\frac{p}{\rho^\gamma})=0
\end{eqnarray}
For sufficiently smooth initial conditions it is irrelevant whether to work with the conservative or the non-conservative formulation.
In case of discontinuities this is not the fact any more.
In a supernova model we have to deal with shock waves, which are characterized by nearly abrupt changes in the physical quantities.
Thus great care has to be taken when forwarding primitive variables numerically.
An example illustrating that the conservative and the non-conservative formulation of a set of PDEs can give different solutions is given by \cite{TOR001} for the case of the shallow water equation.
The eigenvalues of $\frac{\partial f(\mathbf u)}{\partial\mathbf u}$ or $F(\mathbf w)$ are $\lambda_{1,2}=v\pm a$ and $\lambda_{3,4,5}=v=|\mathbf v|$.
All $\lambda_i$ are real, and $F(\mathbf w)$ can be diagonalized.
In one dimension this is $F=P\cdot J\cdot P^{-1}$ with
\begin{eqnarray}
P=\begin{pmatrix}
1 & \frac{1}{a^2} & \frac{1}{a^2} \\
0 & -\frac{1}{a\rho} & -\frac{1}{a\rho} \\
0 & 1 & 1
\end{pmatrix}
\qquad
J=\begin{pmatrix}
v & 0 & 0 \\
0 & v-a & 0 \\
0 & 0 & v+a
\end{pmatrix}
\qquad
P^{-1}=\begin{pmatrix}
1 & 0 & -\frac{1}{a^2} \\
0 & -\frac{a\rho}{2} & \frac{1}{2} \\
0 & \frac{a\rho}{2} & \frac{1}{2}
\end{pmatrix}
\end{eqnarray}
The real eigenvalues with linear independent eigenvectors make the Euler equations a set of hyperbolic PDEs \cite{TRA001}.
\subsubsection{Extending the Euler equations}
The Euler equations are the basis of the model for the expansion of a supernova.
They reflect the flux of a compressible ideal gas with respect to pressure and energy.
However, the solutions of a system of hyperbolic PDEs are generally waves.
If the set of equations is nonlinear, those wave functions are likely to break down at a certain point.
In the case of a supernova expansion shock waves are formed, which result in discontinuities.
The properties of shock waves will be discussed in chapter \ref{SHOCKWAVES}.
In reality physical quantities are never discontinuous.
In the case of the Euler equations, viscosity is responsible for smoothing out the solutions.
Roughly speaking, one has to add some terms to the momentum equation (\ref{EQU002}).
\begin{eqnarray}\label{EQU042}
\partial_t(\rho\mathbf v)+\vec\nabla((\rho\mathbf v)\mathbf v+p\mathbf I)&=&\vec\nabla\cdot\Pi+\mathbf f\label{EQU007}
\end{eqnarray}
Without going into details, $\Pi$ is a tensor of second rank called the {\it viscous stress tensor}, and $\mathbf f$ sums up all kinds of external specific forces acting on the fluid.
The resulting equation is called the {\it Navier--Stokes} equation.
Unfortunately, solving the Navier--Stokes equation can be extremely difficult.
Even making predictions regarding its general behaviour can be very challenging\footnote{It is one of the seven famous millennium problems and a solution is worth US\$\,1,000,000 \cite{BAS001}.}.
In chapter \ref{SECARTIFICALVISCOSITY} we will introduce {\it artificial} viscosity by adding a term proportional to $\rho\frac{\partial v}{\partial r}|\frac{\partial v}{\partial r}|$ to the pressure $p$.
Although this idea was originally suggested by \cite{CIO001} for a Lagrangian approach, it is now a common technique for all kinds of PDEs.
Another extension of the Euler equations is required due to heat conduction.
We will see in chapter \ref{SUPERNOVA} that after a relatively short adiabatic expansion, cooling of the SNR is initiated.
Cooling is a consequence of particle collisions and is therefore proportional to $n^2$ with $n=\frac{\rho}{\mu m_\text{H}}$ being the particle density, $\mu$ the average molecular weight and $m_\text{H}$ the mass of a hydrogen nucleus, i.e.\ a proton.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=\hsize]{SN.fig002.cooling.eps}
\caption{Cooling function: the effect of cooling increases by two orders of magnitude during the cool-down between $10^7\,\mathrm{K}$ and $10^5\,\mathrm{K}$.
This is mainly caused by the presence of carbon and oxygen.}\label{FIG002}
\end{center}
\end{figure}
The cooling function $\Lambda(T)$ describes the radiative loss of a low-density astrophysical plasma.
The data for figure \ref{FIG002} is taken from \cite{SUT001} and nicely illustrates the temperature dependency due to dissociation, recombination and radiation.
Since cooling affects energy, the third Euler equation has to be modified to
\begin{eqnarray}\label{EQU008}
\partial_t\epsilon+\vec\nabla(\epsilon\mathbf v+p\mathbf v)=-n^2\Lambda(T)
\end{eqnarray}
where
\begin{eqnarray}
\frac{d}{dt}\left(\frac{3}{2}nk_BT\right)=\frac{d}{dt}(\rho U)=-n^2\Lambda(T)=-L(T)
\end{eqnarray}
% and $U=\frac{1}{\gamma-1}\frac{p}{\rho}=\frac{c_s^2}{\gamma-1}$ is the specific internal energy for an ideal gas.
and $U=\frac{1}{\gamma-1}\frac{p}{\rho}$ is the specific internal energy for an ideal gas.
$L(T)$ describes the total effect of cooling.
How to implement this for a supernova is discussed for the snowplow phase in chapter \ref{SNOWPLOW}.
\subsubsection{One dimensional Euler equations for different geometries}\label{EULERGEOMETRIES}
In one dimension the Euler equations reduce to a set of 3 equations.
In the following $r$ will be the spatial variable, pointing out the rudimentary symmetry of the supernova model.
Equation (\ref{EQU005}) can be simplified to
\begin{eqnarray}\label{EQU009}
\partial_t\mathbf u+\partial_r f(\mathbf u)=-\frac{\alpha}{r}s(\mathbf u)
\end{eqnarray}
with
\begin{eqnarray}\label{EQU010}
\mathbf u=\begin{pmatrix}\rho\\\rho v\\\epsilon\end{pmatrix}\qquad
f(\mathbf u)=\begin{pmatrix}\rho v\\
\rho v^2+p\\
(\epsilon+p)v
\end{pmatrix}\qquad
s(\mathbf u)=\begin{pmatrix}\rho v\\
\rho v^2\\
(\epsilon+p)v
\end{pmatrix}
\end{eqnarray}
The parameter $\alpha$ characterizes the geometry of the problem, where $\alpha=0$ is the planar case, $\alpha=1$ is the case of cylindrical symmetry and $\alpha=2$ is the case of spherical symmetry.
The term $-\frac{\alpha}{r}s(\mathbf u)$ is called the {\it source term}.
In the planar case ($\alpha=0$) the Euler equations can also be written in a quasi-linear form
\begin{eqnarray}\label{EQU014}
\partial_t\mathbf u+\frac{\partial f(\mathbf u)}{\partial \mathbf u}\partial_x\mathbf u=0,\qquad
\mathbf u=\begin{pmatrix}u_1\\u_2\\u_3\end{pmatrix}=
\begin{pmatrix}\rho\\\rho v\\\epsilon\end{pmatrix},\qquad
f(\mathbf u)=\begin{pmatrix}u_2\\\frac{3-\gamma}{2}\frac{u_2^2}{u_1}+(\gamma-1)u_3\\-\frac{\gamma-1}{2}\frac{u_2^3}{u_1^2}+\gamma\frac{u_2u_3}{u_1}\end{pmatrix}
\end{eqnarray}
Equation (\ref{EQU014}) is the starting point for several numerical schemes.
\subsubsection{Euler equations and supernova}\label{EULERSN}
Consider now the Euler equations in spherical symmetry to simulate the expansion of SNRs.
The different eigenvalues $v$ and $v\pm a$ of the Euler equations correspond to different wave fronts.
Because of the non-linearity of the problem the propagation speeds of these wave fronts are not necessarily proportional to the eigenvalues.
There are only three different possible types of wave fronts \cite{TOR001}
\begin{itemize}
\item shock waves
\item contact waves
\item rarefaction waves
\end{itemize}
The eigenvalue $v$ corresponds to a contact wave, but the eigenvalues $v\pm a$ can principally result in any wave.
In a supernova simulation we have a large initial amount of kinetic energy and therefore $a\ll v$.
This is why both eigenvalues $v\pm a$ are associated with shock waves, the reverse shock for eigenvalue $v-a$ and the shell for $v+a$.
Remember that $a$ denotes as usually the adiabatic speed of sound.
\subsection{Shock waves}\label{SHOCKWAVES}
Shock waves are a very common phenomenon in astrophysics and therefore part of any introduction to this field of knowledge.
I refer to the script of Dorfi \cite{DOR001}.
At a supernova explosion, the velocity of the expansion $v$ is initially 2--3 orders of magnitude larger than the adiabatic speed of sound $a=\sqrt{\frac{\gamma p}{\rho}}$.
This prevents the propagating waves from adapting to the medium.
A thin shock front where the physical quantities change abruptly, is formed.
There are two types of shock waves to distinguish
\begin{itemize}
\item stationary shock waves
\item propagating shock waves
\end{itemize}
An example of stationary shock waves is solar winds.
The shell formed by a supernova explosion on the other hand is a propagating compression shock.
Density, pressure, temperature and energy initially increase, but the propagation velocity of the shock front decreases continuously.
We will formalize and estimate this deceleration in chapter \ref{SUPERNOVA}.
Based on the assumption that no work is done, so that the process is purely adiabatic, the {\it Rankine-Hugoniot} conditions can be derived.
In the rest frame of the shock wave, these can be written \cite{DOR001}\cite{DOR002} as
\begin{eqnarray}
\rho_{\mathrm{L}}\,v_{\mathrm{L}} &=& \rho_{\mathrm{R}}\,v_{\mathrm{R}} \label{EQU011}\\
\rho_{\mathrm{L}}\,v_{\mathrm{L}}^2+p_{\mathrm{L}} &=& \rho_{\mathrm{R}}\,v_{\mathrm{R}}^2+p_{\mathrm{R}} \label{EQU012}\\
\epsilon_{\mathrm{L}}+\frac{p_{\mathrm{L}}}{\rho_{\mathrm{L}}}+\frac{v_{\mathrm{L}}^2}{2} &=& \epsilon_{\mathrm{R}}+\frac{p_{\mathrm{R}}}{\rho_{\mathrm{R}}}+\frac{v_{\mathrm{R}}^2}{2} \label{EQU013}
\end{eqnarray}
Here $\rho_{\mathrm{K}}$, $p_{\mathrm{K}}$, $\epsilon_{\mathrm{K}}$, $v_{\mathrm{K}}$ denote mass density, pressure, total specific internal energy and velocity of the fluid left (L) and right (R) of the shock front.
Similar to the Euler equations, (\ref{EQU011}) embodies conservation of mass, (\ref{EQU012}) conservation of momentum and (\ref{EQU013}) conservation of energy.
Within the framework of this thesis, L will always be the side containing the origin and thus the interior of the supernova and R will be the side of the ISM.
With some calculations \cite{DOR002} one can derive the compression ratio, which is the ratio between the mass density behind the shock and the mass density in front of the shock.
\begin{eqnarray}\label{EQU018}
\frac{\rho_{\mathrm{behind}}}{\rho_{\mathrm{front}}} = \frac{\rho_{\mathrm{shell}}}{\rho_{\mathrm{ISM}}} = \frac{\gamma+1}{\gamma-1}
\end{eqnarray}
For a monoatomic gas $\gamma=\frac{5}{3}$ and the compression ratio is 4:1, for a diatomic gas $\gamma=\frac{7}{5}$ and the compression ratio is 6:1.
In the case of a supernova we are dealing with a hot, monoatomic gas, and the compression ratio is therefore 4:1.
\subsection{Supernova explosion}\label{SUPERNOVA}
Supernovae are extremely fierce events which release extremely large amounts of energy to the interstellar medium (ISM).
They are responsible for a variety of astrophysical processes such as the creation of heavy elements and, more importantly, are essential for stellar evolution.
The explosion itself happens within a fraction of a second.
For studying the development of supernova remnants one has to know that approximately $10^{44}$ Joule (which is as much energy as our sun produces in 7 billion years) are released in a moment, followed by an almost spherically symmetric, self-similar expansion over a few hundred years.
All the released energy is of kinetic nature and the reason for an initial velocity of the ejecta of as much as $5000\,\mathrm{km/s}$ or more.
Basically, there are two types of supernovae.
Type Ia supernovae are caused by carbon deflagration in white dwarfs in binary systems.
There is an upper mass limit for white dwarfs, the {\it Chandrasekhar limit} which is approximately $1.44\,M_\odot$.
By attracting mass from its binary, usually a red giant, the white dwarf gains mass and (counter-intuitively) the radius decreases.
When the Chandrasekhar limit is reached, the radius becomes infinitesimally small, the star collapses and results in a gigantic explosion.
Since the Chandrasekhar limit is in principle always the same, the released energy is nearly the same for all type Ia supernovae.
This fact is utilized for calculating accurate distances of distant objects\footnote{The accelerating expansion of the universe was corroborated with this method \cite{PER001}, Nobel Prize Physics 2011.}.
Supernovae of type II\footnote{There are numerous sub-classifications.} are triggered by the collapse of the core of massive stars in HII-regions\footnote{HII is a an acronym for singly ionized hydrogen, where H stands as usual for hydrogen and one less than the roman number indicates how often the atom is ionized. Thus, HI stands for neutral hydrogen. Similarly, FeVI stands for fivefold ionized iron.}, the origin of most stars.
Depending on their mass they leave a neutron star or a black hole behind.
A typical number for mass ejected by such a supernova is $3\,M_\odot$ upwards.
Massive main sequence stars (OB-associations) have a short lifetime not exceeding a few million years.
A type II supernova is therefore likely to occur close to the birthplace of the star, thus in a HII region.
In the present simulation we will focus our attention on type II supernovae.
We will assume the hot environment of an HII region and investigate the interaction with the surrounding medium.
In chapter \ref{ISM}, we will discuss the properties of this neighborhood.
The expansion of the supernova remnants will be explained in chapter \ref{SNR}.
\subsection{Interstellar medium, HII regions}\label{ISM}
The physical conditions in the vicinity of the supernova play a major role.
It is therefore absolutely necessary to know the nature and the physical processes of the ISM in a HII region.
A typical HII region has an extension of 100--1000 pc and is the birthplace of most stars.
It contains plenty of young, hot, blue stars.
While most of a HII region consists of hot, ionized hydrogen, every now and then local molecular condensations develop, cool down¸ collapse and form new stars.
The ISM is, in simplest approximation, a two-phase model \cite{FIE001} with a cold phase at $T=80\,\mathrm{K}$ and a warm phase at $T=8000\,\mathrm{K}$.
The reasons for this are found in the equilibrium between the heating and cooling processes and the second law of thermodynamics \cite{DOR002}.
A third phase consists of the slowly cooling SNRs, a very hot, thin gas at temperatures $>10^{6}\,\mathrm{K}$.
This phase is not in equilibrium with either of the other two phases and is nourished by constant new supernova explosions.
During the life span of an HII region of a few hundred Myrs, thousands of stars are born.
Of these only the extremely massive OB-associations are of interest to us.
The average lifetime of an OB-association strongly depends on its mass and is between $0.5\,\mathrm{Myrs}$ for $40\,M_\odot$ O-stars and $100\,\mathrm{Myrs}$ for $5\,M_\odot$ B-stars.
The very massive O-stars are rare and therefore, for an OB-association one can reasonably assume a life span of approximately $\tau=50\,\mathrm{Myrs}$.
Based on observations it is estimated that during a period $\tau$ about N=30--50 supernova explosions \cite{DOR002} occur.
This results in an average power $L_\mathrm{SN}$ produced by supernovae of
\begin{eqnarray}
L_\mathrm{SN}=\frac{N\cdot 10^{44}\mathrm{W}}{\tau}\approx 2\cdot 10^{30}\mathrm{W}
\end{eqnarray}
This is four magnitudes higher then the luminosity of our sun $L_\odot=3.8\cdot 10^{26}\,\mathrm{W}$ and one of the reasons for the warm phase at temperatures around $8000\,\mathrm{K}$.
The temperature is an important parameter in the model because it is proportional to the internal pressure $p$ and the energy density $\rho U$.
\begin{eqnarray}\label{EQU015}
p=n k_B T=\frac{\rho}{\mu m_\mathrm{H}}k_B T
\end{eqnarray}
The other important parameter in (\ref{EQU015}) is the number density $n$ of the environment.
More significantly, in the hot phase, the number density is according to \cite{MCK001} somewhere between $n=0.25\,/\mathrm{cm}^3$ and $n=0.37\,/\mathrm{cm}^3$.
In order to study the interaction between the ISM and a SNR, our later model will be based on a particle density of $n=1/\mathrm{cm}^3$.
The cold phase precedes the gravitational collapse that happens prior to star birth.
In the cold phase, at $T=80\,\mathrm{K}$, the particle density according to \cite{MCK001} is typically $n=42\,/\mathrm{cm}^3$.
We will see later that this density is almost one order of magnitude higher than the density of the SNR after the initial free expansion of a few hundred years.
For these reasons we will not consider the cold phase in our models.
Finally $\mu$ is the mean molecular weight.
For a fully ionized gas containing only hydrogen and helium this is
\begin{eqnarray}
\mu=\frac{4-X}{3(X+1)}=\frac{Y+3}{3(2-Y)}
\end{eqnarray}
where $X$ is the mass fraction of hydrogen and $Y=1-X$ is the mass fraction of helium respectively.
Using $X=0.75\,/\,Y=0.25$ this gives $\mu=\frac{13}{21}\approx 0.62$.
\subsection{Supernova remnants (SNR)}\label{SNR}
The evolution of SNRs and their interaction with the ISM are essential for our supernova model.
In the following pages we derive a construct for the theoretical description of this process.
These findings are used to control and ultimately to fine-tune the numerical simulations.
There are four major stages during the expansion of a supernova:
\begin{itemize}
\item initial phase: free expansion
\item Sedov-Taylor phase: adiabatic expansion
\item snowplow phase: radiative heat conduction of the shell
\item final phase: cooling of the interior of the supernova and mixture with the ISM
\end{itemize}
To make calculations easier we introduce a set of dimensionless quantities.
\begin{align}
& M^\ast = \frac{M}{M_\odot} &
& R^\ast = \frac{R}{\mathrm{pc}} &
& t^\ast = \frac{t}{10^4\,\mathrm{yrs}} & \\
& n^\ast = \frac{n}{10^6\,/\mathrm{m}^3} &
& \rho^\ast = 0.0247\mu\, n^\ast &
& E^\ast_\mathrm{SN} = \frac{E_\mathrm{SN}}{10^{44}\,\mathrm{J}} = 5258.398^{-1}\,\frac{\mathrm{M}^\ast{\mathrm{R}^\ast}^2}{{\mathrm{t}^\ast}^2} \\
% & \rho^\ast = \frac{M^\ast}{{R^\ast}^3} &
% & p^\ast = \frac{M^\ast}{R^\ast{t^\ast}^2} &
% & E^\ast_{44} = 5258.4\,\frac{E_\mathrm{SN}}{E_{44}} &
\end{align}
Here $E_\mathrm{SN}$ is the total released energy by the supernova explosion.
\paragraph{Free expansion:}
The {\it sweep-up radius} $R_\mathrm{s}$ is defined as the radius, where the displaced mass of the ISM becomes equal to the mass of the supernova ejecta $M$.
\begin{eqnarray}
M &\overset{!}{=}& \frac{4\pi}{3}\rho_\mathrm{ISM}R_\mathrm{s}^3 = \frac{4\pi}{3}\frac{n}{\mu m_\mathrm{H}} R_\mathrm{s}^3\\
\Longrightarrow\qquad R^\ast_\mathrm{s} &=& 2.130\times\sqrt[3]{\frac{M^\ast}{\mu n^\ast}}\label{EQU016}
\end{eqnarray}
The total energy $E_\mathrm{SN}$ of the supernova is converted into kinetic energy and the initial velocity $v_\mathrm{SN}$ of the ejecta is therefore approximately
\begin{eqnarray}\label{EQU017}
v_\mathrm{SN}=\sqrt{\frac{2E_\mathrm{SN}}{M}}
\end{eqnarray}
During the free expansion, the position $R(t)=v_\mathrm{SN}t$ of the shock wave increases almost linearly.
This phase lasts until the sweep-up time $t_\mathrm{s}=R_\mathrm{s}/v_\mathrm{SN}$, when the shock wave reaches the sweep-up radius $R_\mathrm{s}$.
Using (\ref{EQU016}) and (\ref{EQU017}) this gives
\begin{eqnarray}
t^\ast_\mathrm{s}=0.0208\,{M^\ast}^{\frac{5}{6}}\,({\mu n^\ast})^{-\frac{1}{3}}\,{E^\ast_\mathrm{SN}}^{-\frac{1}{2}}
\end{eqnarray}
with typical values of a few hundred years.
Our supernova model starts at some point during the free expansion.
From then on the behaviour of the SNR is characterized by the solution of the Euler equations.
Because of the high velocity of the ejecta, a shock wave is formed immediately.
The smallest eigenvalue of the Euler equations is responsible for a reverse shock, which starts travelling inwards and is a well known phenomenon.
Choosing proper initial conditions is not a trivial task.
We will pursue two different ways.
The first one is derived from a rough estimate of the density profile and will be explained below.
The second one is more sophisticated and uses the exact blast wave solution given by \cite{SED002}.
While the first approach provides some freedom in the parameters, the second one is restrictive and will be outlined separately in chapter \ref{SEDOVSOLUTION}.
For the simpler approach we assume a self-similar expansion in a linear velocity field for the region between the origin and the shell
\begin{eqnarray}
v=\frac{r}{t}=\frac{\partial r}{\partial t}
\end{eqnarray}
To satisfy the first Euler equation in spherical coordinates
\begin{eqnarray}
\partial_t\rho+\frac{1}{r^2}\partial_r(r^2v\rho)=0
\end{eqnarray}
one uses the ansatz \cite{DWA002}
\begin{eqnarray}
\rho(r,t)=\rho(v,t)=\frac{A}{t^3}e^{-\frac{v}{v_0}}
\end{eqnarray}
where the constants $A$ and $v_0$ are determined by mass and energy conservation.
\begin{align}
& 4\pi\int_0^\infty \rho r^2\,dr = M &
& 4\pi\int_0^\infty \frac{\rho v^2}{2}r^2 dr = E & \\
& A = \frac{M}{8\pi}v_0^{-3} &
& v_0 = \sqrt{\frac{E}{6M}} &
\end{align}
This leads to an initial density profile
\begin{eqnarray}\label{EQU040}
\rho(r,t)=\frac{6^\frac{3}{2}}{8\pi}\,M^\frac{5}{2}\,E^{-\frac{3}{2}}\,t^{-3}\, e^{-\sqrt{\frac{6M}{E}}\frac{r}{t}}
\end{eqnarray}
At this initial point, we assume the shock wave already to be at position $r_0$.
The position $r_0$ should be far enough from the origin to avoid numerical problems, but close enough to describe processes such as the reverse shock.
The starting time $t_0$ of the simulation is consequently at
\begin{eqnarray}
t_0=\frac{r_0}{v_\mathrm{SN}}
\end{eqnarray}
A reasonable number for $r_0$ could be $0.5\,\mathrm{pc}$ and $t_0$ could be roughly 80 years after the blast.
Another essential initial ingredient is the pressure distribution.
While the pressure $p$ of the ISM in front of the shock wave is explicitly given by the density $\rho_\mathrm{front}=\rho_\mathrm{ISM}$ and the temperature $T_\mathrm{ISM}$ of the ISM
\begin{eqnarray}
p=\frac{\rho_\mathrm{ISM}k_B T_\mathrm{ISM}}{\mu m_\mathrm{H}}
\end{eqnarray}
the situation behind the shock wave needs some further considerations.
One way is to use the Rankine-Hugoniot equations (\ref{EQU011}) and (\ref{EQU012}) and the compression ration (\ref{EQU018}).
We neglect the pressure $p_\mathrm{ISM}=p_{\mathrm{R}} \ll p_{\mathrm{L}}$ in front of the shock, use $v_{\mathrm{R}}=v=\frac{dr}{dt}$ for the velocity of the shock wave and set $\rho_{\mathrm{R}}=\rho_\mathrm{ISM}$ for the density of the ISM.
Taking all this into account one gets
\begin{eqnarray}\label{EQU023}
p_{\mathrm{L}}=p_\mathrm{behind}=\frac{2\rho_\mathrm{ISM}(\frac{dr}{dt})^2}{\gamma+1}
\end{eqnarray}
for the pressure behind the shock.
Unfortunately this is only the pressure close behind the shock wave and does not give any further information about the full pressure profile towards the origin.
For the approach described here it is sufficient to use a piecewise constant pressure profile
\begin{eqnarray}
p_\mathrm{init} &=& \begin{cases}
p_{\mathrm{L}}=p_\mathrm{behind} & \qquad r\le r_0 \\
p_{\mathrm{R}}=p_\mathrm{ISM} & \qquad r> r_0 \\
\end{cases}
\end{eqnarray}
% For the rough model we can set for $0\le r\le r_0$ the pressure $p_\mathrm{behind}$ constant according to (\ref{EQU023}).
The simulations will show that a reasonable pressure profile develops immediately.
However, in chapter \ref{SEDOVSOLUTION} an exact blast wave solution for more accurate simulations will be given.
\paragraph{Sedov-Taylor phase:}
As soon as the mass of the supernova becomes negligible compared to the displaced mass of the ISM, a mainly pressure-driven adiabatic blast wave phase starts.
We will again start with a simple approximation and refer to the exact self-similar solution as given in \cite{TAY001} in chapter \ref{SEDOVSOLUTION}.
Luckily this simpler approach gives an almost correct result.
Assume that all collected matter from the ISM at time $t$ is located in a very thin shell at a distance $r(t)$ from the origin.
This assumption makes sense as can be seen by the following argument:
let $\delta r$ be the radial size of the shell, $\rho_\mathrm{behind}$ the density immediately behind the shock wave and $\rho_\mathrm{front}=\rho_\mathrm{ISM}$ the density of the ISM in front of the shock wave.
Using (\ref{EQU018}) one estimates the spreading of the shell for an ideal monoatomic gas ($\gamma=\frac{5}{3}$) by
\begin{eqnarray}
4\pi r^2\delta r\rho_\mathrm{behind}&=&\frac{4\pi}{3}r^3\,\rho_\mathrm{front}\\
\delta r = \frac{\rho_\mathrm{front}}{\rho_\mathrm{behind}}\cdot\frac{r}{3} = \frac{\gamma-1}{\gamma+1}\cdot\frac{r}{3} &=& \frac{r}{12}
\end{eqnarray}
So the extent of the shock wave is relatively small compared to the expansion radius.
In our rough approximation we consider $\delta r$ to be infinitesimally small.
The pressure in the shell must compensate the expanding matter and therefore
\begin{eqnarray}\label{EQU019}
4\pi r(t)^2 p(r(t)) &=& \frac{d}{dt}\left(M(t)v(t)\right)=\frac{d}{dt}\left(\frac{4\pi}{3}r(t)^3 \rho_\mathrm{ISM} \frac{dr(t)}{dt}\right)
\end{eqnarray}
Finally we substitute $p(t)$
\begin{eqnarray}
(\gamma-1)E_\mathrm{SN} &=& p(t)V(t) = \frac{4\pi}{3}r(t)^3p(t)
\end{eqnarray}
and use the ansatz
\begin{eqnarray}
r(t)=A\cdot t^\alpha
\end{eqnarray}
which leads to
\begin{eqnarray}\label{EQU020}
r(t)=\left(\frac{75(\gamma-1)E_\mathrm{SN}}{8\pi\rho_\mathrm{front}}\right)^\frac{1}{5}\cdot t^\frac{2}{5}
\end{eqnarray}
Rewriting (\ref{EQU020}) in dimensionless form gives
\begin{eqnarray}\label{EQU021}
{r^\ast} &=& 14.473\cdot\sqrt[5]{\frac{(\gamma-1)E^\ast {t^\ast}^2}{\mu n^\ast}}
\end{eqnarray}
Equation (\ref{EQU021}) will be one essential criterion to check and fine-tune the supernova model.
The exact self-similar solution in chapter \ref{SEDOVSOLUTION} using $\gamma=\frac{5}{3}$ and $\mu=\frac{13}{21}$ is bigger by just about $0.3\%$.
\begin{eqnarray}\label{EQU022}
% {r^\ast}-r^\ast_\mathrm{s} &=& 14.7\cdot\sqrt[5]{\frac{E^\ast (t^\ast-t^\ast_\mathrm{s})^2}{n^\ast}}
{r^\ast} &=& 14.7421\cdot\sqrt[5]{\frac{E^\ast {t^\ast}^2}{n^\ast}}
\end{eqnarray}
% for adjusting the model, also taking sweep-up radius $r^\ast_\mathrm{s}$ and sweep-up time $t^\ast_\mathrm{s}$ into account.
Obviously, this formula makes sense only for $t^\ast>t^\ast_\mathrm{s}$.
There is one more essential issue to be addressed.
After the linear expansion, the further outspread of the SNR obeys the Euler equations.
In chapter \ref{EULER} we have learned that there are three eigensolutions corresponding to the three eigenvalues $v$ and $v\pm a$.
The eigensolution with eigenvalue $v-a$ is responsible for a reverse shock, which initially runs inwards and heats up the interior of the SNR.
This reverse shock can be observed in the simulations in chapter \ref{SEC005}.
However, some precautions have to be taken since this reverse shock can travel back to the origin, where the coordinate singularity causes numerical problems.
% The temporal development of the temperature $T(t)$ of the shell during the Sedov-Taylor phase can be found e.g.\ in \cite{DOR002}.
\paragraph{Snowplow-phase:}\label{SNOWPLOW}
In the beginning the temperatures in the shell are extremely high.
While the interior of the SNRs continues to expand adiabatically, the shell starts to cool-down.
Figure \ref{FIG002} shows nicely how the efficiency of cooling increases rapidly at temperatures below $10^6\,\mathrm{K}$.
At these temperatures the cooling follows approximately \cite{FIE001}
\begin{eqnarray}\label{EQU024}
\Lambda(T)\propto\sqrt{\frac{1}{T}}
\end{eqnarray}
Before the shell enters the phase of radiative cooling, the temperature in it can be estimated by the ideal gas equation, (\ref{EQU023}) and (\ref{EQU020}).
\begin{eqnarray}
p &=& \frac{n}{V}k_B T = \frac{\rho_\mathrm{behind}}{m_\mathrm{H}\mu}k_B T=\frac{\frac{\gamma+1}{\gamma-1}\rho_\mathrm{ISM}}{m_\mathrm{H}\mu}k_B T=
\frac{2\rho_\mathrm{ISM}v_\mathrm{shell}^2}{\gamma+1}\\
k_B T &=& \frac{2(\gamma-1)m_\mathrm{H}\mu}{(\gamma+1)^2}\cdot\left(\frac{75(\gamma-1)E_\mathrm{SN}}{8\pi\rho_\mathrm{ISM}}\right)^\frac{2}{5}\cdot\frac{4}{25}t^{-\frac{6}{5}}\\
T(t^\ast) &=& 4640496\cdot\left(\frac{E^\ast_\mathrm{SN}}{n^\ast}\right)^{\frac{2}{5}}\cdot {t^\ast}^{-\frac{6}{5}}\,\mathrm{K}
\end{eqnarray}
Here $v_\mathrm{shell}=dr/dt$ is the propagation velocity of the shock front and we have used $\gamma=\frac{5}{3}$ and $\mu=\frac{13}{21}$ as before.
Therefore, after about 30000--40000 years the temperature of the shell has fallen below $10^6\,\mathrm{K}$ and cooling is well under way.
In \cite{DOR002} there is a derivation for the time $t_\mathrm{cool}$ where the age of the supernova $t_\mathrm{dyn}$ is about the same as the average cooling time $t_c$, which is approximately the internal energy $\epsilon$ divided by the average cooling rate $L$.
\begin{eqnarray}
t_\mathrm{dyn}\approx t_\mathrm{c}\approx\frac{\epsilon}{L}=\frac{\epsilon}{n^2\Lambda(T)}
\end{eqnarray}
The result, however, is dependent on the
% inaccurate
proportionality constant in (\ref{EQU024}) and I shall therefore confine myself to quoting the result
\begin{eqnarray}
t^\ast_\mathrm{cool} &\approx& 2.13\cdot{n^\ast}^{-\frac{4}{7}}{E_\mathrm{SN}^\ast}^{\frac{3}{14}}\\
r^\ast_\mathrm{shell}(t^\ast_\mathrm{cool}) &\approx& 20\cdot{n^\ast}^{-\frac{3}{7}}{E_\mathrm{SN}^\ast}^{\frac{2}{7}}
\end{eqnarray}
The point in time $t^\ast_\mathrm{cool}$ is to an extent of interest to us because here a slowing of the propagation velocity takes place, which is observed in the model simulations.
As mentioned above, the interior of the SNRs continues to expand adiabatically for quite some time, i.e.
\begin{eqnarray}
pV^\gamma\equiv\mathrm{const}
\end{eqnarray}
The ram--pressure $p_\mathrm{ram}=p$ is proportional to the square of the propagation velocity $v_\mathrm{shell}$ and, therefore, for $\gamma=\frac{5}{3}$
\begin{eqnarray}
\dot r^2 r^5 &\equiv& \mathrm{const}\\
r(t) &=& r_\mathrm{shell}(t_\mathrm{cool})\left(\frac{t}{t_\mathrm{cool}}\right)^{\frac{2}{7}}
\end{eqnarray}
For the simulations the cooling function needs to be incorporated into the Euler equations.
This is done by adding the cooling function to the energy equation
\begin{eqnarray}
\partial_t\epsilon+\vec\nabla(\epsilon\mathbf v+p\mathbf v)&=&-\left(\frac{\rho}{\mu m_\text{H}}\right)^2\Lambda(T(\rho,p))\label{EQU025}
\end{eqnarray}
with
\begin{eqnarray}
T(\rho,p) &=& \frac{\mu m_\text{H}}{k_B}\cdot\frac{p}{\rho}
\end{eqnarray}
\paragraph{Final stage:}
The expansion of the SNR will last until its expansion velocity becomes comparable with the typical velocity of the ISM.
This will be after a few Myrs.
A supernova explosion is usually not an isolated event, but happens in a HII region, where thousands of stars are born within a short period.
The shock fronts of adjacent supernovae will meet and the assumption of spherical symmetry will get lost.
Through the collision of SNRs, regions of higher density will develop and alternate with regions of lower density.
These are the conditions required for further star generation.
The consequence for our model is that at this point we have to give up spherical symmetry in the differential schemes.
\newpage
\subsubsection{Summary of parameters}
The following table gives typical values for the parameters.
\begin{center}
\begin{tabular}{r|c|l}
parameter & typical values & description \\
\hline
$E_\mathrm{SN}$ & $10^{44}\,\mathrm{J}$ & energy released by SN \\
$M$ & $3\,M_\odot$ & total mass of ejecta \\
$v_\mathrm{SN}$ & $5789\,\mathrm{km/s}$ & average initial velocity of ejecta \\
$\gamma$ & $\frac{5}{3}$ & adiabatic index for monoatomic gas \\
$X$ & $0.75$ & mass fraction hydrogen in ISM \\
$Y$ & $0.25$ & mass fraction helium in ISM \\
$\mu$ & $\frac{13}{21}\approx 0.62$ & mean molecular weight \\
$n$ & $1/\mathrm{cm}^3$ & particle number density of ISM \\
$\rho_\mathrm{ISM}$ & $1.035\cdot 10^{-21}\,\mathrm{kg/m^3}$ & mass density of ISM \\
$T_\mathrm{ISM}$ & $10000\,\mathrm{K}$ & temperature of ISM \\
$t_\mathrm{s}$ & $610\,\mathrm{years}$ & sweep-up time \\
$R_\mathrm{s}$ & $3.6\,\mathrm{pc}$ & sweep-up radius \\
$t_\mathrm{cool}$ & $21300\,\mathrm{years}$ & time when cooling becomes relevant \\
$R_\mathrm{cool}$ & $20\,\mathrm{pc}$ & position of shell when cooling becomes relevant \\
$T_\mathrm{cool}$ & $2\cdot 10^6\,\mathrm{K}$ & temperature of shell when cooling becomes relevant \\
\end{tabular}
\end{center}
\newpage
\section{Mathematical background}\label{SEC003}
\begin{abstract}
{\it
This chapter covers mathematical and numerical issues.
It starts with general notes on numerically solving PDEs, such as the monotonicity preserving conditions given by Godunov.
The Riemann problem of solving the Euler equations with one initial discontinuity and otherwise constant initial values will be introduced.
The exact solution of the Riemann problem is not only required for a number of differential schemes, but is also a useful tool for verifying numerical solutions.
There is a long list of differential schemes introduced and advantages and disadvantages are outlined.
To get plausible initial conditions shortly after the blast, two self-similar solutions for a punctual explosion as proposed by Taylor and Sedov are given.
}
\end{abstract}
\subsection{Solving PDEs numerically}\label{PDENUMERICALLY}
% In der numerischen Mathematik tauchen Riemann-Probleme in natürlicher Weise in Finite-Volumen-Verfahren zur Lösung von Erhaltungsgleichungen auf.
Solving systems of nonlinear partial differential equations numerically is not easy!
Even comparatively simple problems, such as the inviscid Burgers' equation $\partial_t\mathbf{u}+\mathbf{u}\partial_x\mathbf{u}=0$, immediately reveal the full range of problems that might occur.
It is naive to presume that accuracy can be gained by finer discretization.
A better fit due to smaller distances between the nodes will be more than compensated by the need for shorter time steps.
As already mentioned in chapter \ref{EULERSUB} even the formulation of the problem can be essential.
Using a set of conservative rather than non-conservative quantities gives different solutions in the case of non-continuous initial values (see counterexample in \cite{TOR001}).
The original objective of this thesis was to find out which scheme should be used for which problem.
Before possible approaches are introduced in chapter \ref{SCHEMES}, some prerequisites are required.
All numerical PDEs discussed here have the common form
\begin{eqnarray}
\mathbf u_i^{(n+1)} = \sum_{j=-j_{\mathrm{L}}}^{j_{\mathrm{R}}}b_j\mathbf u_{i+j}^{(n)}
\end{eqnarray}
where the index $i$ is the spatial component and $(n)$ indicates the $n$-th iteration in time.
It is obvious that the $b_j$ depend on the fineness of the spatial grid $\Delta x$ as well as on the time step $\Delta t$.
We talk about a monotone scheme if all the coefficients $b_j\ge0$.
% In 1959, Godunov proved the following fundamental theorem\footnote{A comprehensive proof is given in Wikipedia.}:
In 1959, Godunov proved the following fundamental theorem \cite{TOR001}:
\begin{theorem}
Monotone schemes are of first order only.
\end{theorem}
This is bad news since non-monotone schemes tend to produce oscillations over time.
So the natural problem is to diminish these oscillations by calculating the $b_j$ in an adequate way.
In the following we assume that the system of PDEs to be solved has the one-dimensional form
\begin{eqnarray}
\partial_t\mathbf u(t,x)+\partial_x f(\mathbf u(t,x))=g(\mathbf u(t,x))
\end{eqnarray}
The Euler equations are precisely of this form.
The term $g(\mathbf u)$ is sometimes called {\it source term} and can either originate from the geometry of the problem or from additional terms such as heat conduction or viscosity.
In any case, $g(\mathbf u)$ only depends on $\mathbf u$ but not on its derivatives.
With this, the problem can be split into a homogeneous and an inhomogeneous part, which have to be solved separately.
For the homogeneous part $\partial_t\mathbf u+\partial_x f(\mathbf u)=0$ numerical schemes are used, which in principle follow the pattern
\begin{eqnarray}
\mathbf u_i^{(n+1)}=\mathbf u^{(n)}_i-\frac{\Delta t}{\Delta x}\left(\mathbf f_{i+\frac{1}{2}}-\mathbf f_{i-\frac{1}{2}}\right)
\end{eqnarray}
The $\mathbf f$ are fluxes.
Half-indices indicate that these fluxes are calculated between the grid points\footnote{Actually the vectorial fluxes are located on integer coordinates while the scalar values are positioned somewhere between.
For finite volume methods this is usually at the mass center of the volume.}.
Several decades of research has given us suitable techniques to estimate these fluxes appropriately.
% It is the result of decades of research that there are suitable techniques to estimate these fluxes appropriately.
Once the homogeneous part is done and $\mathbf u$ is forwarded by one time step $\Delta t$, an ordinary differential equation is left
\begin{eqnarray}
\frac{d\mathbf u^{(n+1)}}{dt}=g(\mathbf u^{(n+1)})
\end{eqnarray}
This ODE can be forwarded by simple Runge--Kutta methods or other suitable algorithms.
Not every consistent differential scheme is numerically stable.
In 1947 the Hungarian mathematician Von Neumann developed a framework for linear problems, known as {\it Von Neumann stability analysis}, to decide, whether a scheme is stable or not.
The basic idea is, to decompose numerical errors into Fourier modes.
If the (absolute value) of an eigenvalue of one of these modes exceeds unity, the scheme will become unstable after a finite time.
It is important to understand (\smiley) that this analysis deals with the effects of numerical errors rather than errors caused by discretization.
For the Euler equations, stability analysis for the different schemes is substance for extensive further investigations.
Here an easier way is adopted.
The stability of the various differential schemes is explored by some simple test cases described in chapter \ref{SEC004}.
The useful schemes are then applied on the Euler equations and the unstable ones are discarded.
Basically, there are four fundamental properties of numerical schemes \cite{TOR001}
\begin{itemize}
\item consistency
\item stability
\item accuracy
\item convergence
\end{itemize}
{\it Consistency} means that effectively the correct system of PDEs is solved, i.e.\ the differential scheme approximates the original set of equations.
That this is not necessarily always the case is shown by the counterexample for the Du Fort-Frankel-scheme applied to the heat equation given in \cite{NEU003}(p64).
The {\it stability} of a scheme is established, if its fundamental convergence behavior is not affected by rounding errors, independent from the numerical quality.
On the other hand, {\it accuracy} describes the difference between the exact and the numerical solution.
It is defined \cite{TOR001}(p417) as the leading term in the local truncation error.
{\it Convergence} is hard to prove.
At least for linear systems the {\it Lax equivalence theorem} states for well-posed problems
\begin{theorem}
In linear systems consistency and stability is equivalent with convergence.
% TORO, p418
\end{theorem}
Since consistency is easy to see, in the linear case one only has to ensure that a given scheme is stable to converge.
Many differential schemes require solving a countless number of isolated Riemann problems.
Therefore we will give a short discussion of the Riemann problem in chapter \ref{RIEMANNPROBLEM} before we advance to numerical schemes in chapter \ref{SCHEMES}.
\subsection{Riemann problem}\label{RIEMANNPROBLEM}
A Riemann problem is an initial value problem with constant initial values, except for a hypersurface where there is a discontinuity.
Riemann problems are ubiquitous and describe phenomena such as shock waves or rarefaction waves.
It is surprising at first sight that for the one-dimensional Euler equations (which are hyperbolic and nonlinear), exact solutions of the Riemann problem can be constructed.
However, in practice, solutions are approximated by so-called Riemann solvers.
A large number of such numerical Riemann solvers are given in \cite{TOR001}.
The only Riemann solver that will be touched upon is the one by Roe \cite{ROE001}, which is quite common in hydrodynamic simulations and finds its application e.g.\ in the simulations of \cite{SHA001} or \cite{NIS001}.
In a very simple approximation even the supernova model can be assumed to be a Riemann problem.
The initial point of the shell $r_0$ marks the discontinuity, and the total mass of the ejecta is uniformly distributed within the shell.
\begin{align}
\rho_{\mathrm{L}}=\frac{3M}{4\pi r_0^3} &\qquad \rho_{\mathrm{R}}=\rho_{\mathrm{ISM}} &\qquad \text{initial density} \\
v_{\mathrm{L}}=\sqrt{\frac{2E_\mathrm{SN}}{M}} &\qquad v_{\mathrm{R}}=0 &\qquad \text{initial velocity} \\
p_{\mathrm{L}}=\frac{2\rho_{\mathrm{ISM}}v_{\mathrm{L}}^2}{\gamma+1} &\qquad p_{\mathrm{R}}=n_{\mathrm{ISM}}k_B T_{\mathrm{ISM}} &\qquad \text{initial pressure} \\
\end{align}
Such initial conditions are only a very rough approximation.
However, after a couple of time steps the general pattern of the solution becomes similar to more accurate models.
Here we have to presuppose spherical symmetry, and unfortunately there is no exact solution of the Riemann problem in this case.
In this thesis, the exact solution of the Riemann problem is exploited in two complementary ways:
\begin{enumerate}
\item As a sub-tool for numerical schemes:\\
due to discretization, the physical quantities during a simulation can be considered as piecewise constant within a range with discontinuities between the finite volumes.
Many numerical schemes therefore require successively the solutions of a large number of Riemann problems.
For every time step and spatial point a Riemann problem is solved and the individual results are compiled to an overall solution.
\item As an instrument to verify numerical results:\\
the exact solution of a Riemann problem is one way to verify numerical schemes.
We will later use the exact solution of a plane, one-dimensional Riemann problem and compare it to different simulation results.
\end{enumerate}
\subsubsection{Exact solution of the Riemann problem for the one-dimensional plane Euler equations}
The Riemann problem is stated as the solution of the Euler equations with constant initial conditions except one discontinuity.
We will consider the planar one-dimensional case\footnote{The spherically symmetric case involves a source term and can not be treated so easily.}, thus
\begin{eqnarray}\label{EQU041}
\partial_t\mathbf u+\vec\nabla f(\mathbf u)=0,\qquad
\mathbf u=\begin{pmatrix}\rho\\\rho v\\\epsilon\end{pmatrix},\qquad
f(\mathbf u)=\begin{pmatrix}\rho v\\
\rho v^2+p\\
(\epsilon+p) v
\end{pmatrix}
\end{eqnarray}
with initial conditions
\begin{eqnarray}
\mathbf u(x,0)=
\begin{cases}
\begin{pmatrix}\rho_{\mathrm{L}}\\\rho_{\mathrm{L}} v_{\mathrm{L}}\\\frac{1}{2}\rho_{\mathrm{L}} v_{\mathrm{L}}^2+\frac{p_{\mathrm{L}}}{\gamma-1}\end{pmatrix}\qquad(xx_0)
\end{cases}
\end{eqnarray}
and EOS
\begin{eqnarray}
p=(\gamma-1)(\epsilon-\frac{1}{2}\rho v^2)
\end{eqnarray}
Note that we are solving the conservative form of the Euler equations for $\mathbf u=(\rho,\rho v,\epsilon)$ using the non-conservative variables $\mathbf w=(\rho,v,p)$.
As already outlined in chapter \ref{EULERSUB} there are three eigenvalues $\lambda_1=v-a$, $\lambda_2=v$ and $\lambda_3=v+a$ each corresponding to a nonlinear wave front.
These waves split the $(x-t)$-plane into the 4 sections \rm{I}, \rm{III}, \rm{IV} and \rm{VI} as illustrated in figure \ref{FIG001}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=\textwidth]{SN.fig034.eps}
\caption{The wave fronts corresponding to the three eigenvalues split the $(x-t)$-plane into up to six domains.
At the border between \rm{III} and \rm{IV} the density $\rho$ takes a step, but the velocity $v$ and the pressure $p$ remain constant.}\label{FIG001}
\end{center}
\end{figure}
The wave front between \rm{III} and \rm{IV} is always a contact discontinuity.
There is a sudden step in the density $\rho$, but velocity $v$ and pressure $p$ remain constant.
The other two waves associated with the eigenvalues $v\pm a$ can either be shock waves or rarefactions.
In figure \ref{FIG001}, the left wave is a rarefaction and the right wave is a shock wave, but it is also possible that both are shock waves or both are rarefactions.
In the case of a rarefaction, the $(x-t)$-plane is further divided, taking into account the expansion of the rarefaction.
Here \rm{II} is the spread of the rarefaction, while \rm{V} is the infinitesimally small width of the shock wave and can be ignored.
Concerning supernovae, shock waves are far more important than rarefactions.
Moreover, shock wave solutions are slightly easier since one has to treat one case less, namely the expansion of the rarefaction.
The following chapter \ref{SWSOL} outlines the exact approach for both cases and gives the full solution for shock waves explicitly.
\subsubsection{Shock wave solution}\label{SWSOL}
There are two conditions for the occurrence of shock waves \cite{TOR001}, which are
\begin{eqnarray}
\lambda(\mathbf u_{\mathrm{R}})<&\sigma&<\lambda(\mathbf u_{\mathrm{L}})\label{EQU036}\\
f(\mathbf u_{\mathrm{R}})-f(\mathbf u_{\mathrm{L}})&=&(\mathbf u_{\mathrm{R}}-\mathbf u_{\mathrm{L}})\sigma\label{EQU037}
\end{eqnarray}
Here $\lambda$ denotes the eigenvalue, $\sigma$ the propagation speed of the wave and the vector $f(\mathbf u)$ is given in equation (\ref{EQU041}).
The first condition is the Lax entropy condition which manifests the compressive character of the wave while the second condition gives the Rankine-Hugoniot equations when transforming into the rest frame of the wave, thus $\sigma=0$ and $v'_{\mathrm{L}}=v_{\mathrm{L}}-\sigma$, $v'_{\mathrm{R}}=v_{\mathrm{R}}-\sigma$.
\begin{eqnarray}
\rho_{\mathrm{L}} v'_{\mathrm{L}} &=& \rho_{\mathrm{R}} v'_{\mathrm{R}}\label{EQU033}\\
\rho_{\mathrm{L}} (v'_{\mathrm{L}})^2+p_{\mathrm{L}} &=& \rho_{\mathrm{R}} (v'_{\mathrm{R}})^2+p_{\mathrm{R}}\label{EQU034}\\
\frac{\gamma}{\gamma-1}\frac{p_{\mathrm{L}}}{\rho_{\mathrm{L}}}+\frac{1}{2}(v'_{\mathrm{L}})^2 &=& \frac{\gamma}{\gamma-1}\frac{p_{\mathrm{R}}}{\rho_{\mathrm{R}}}+\frac{1}{2}(v'_{\mathrm{R}})^2\label{EQU035}
\end{eqnarray}
In equations (\ref{EQU036})-(\ref{EQU035}) the indices $\mathrm{L}$ and $\mathrm{R}$ stand for left and right, respectively, of a particular shock wave.
Let us now consider the Riemann problem with a two shock wave solution.
The indices $\mathrm{L}$ and $\mathrm{R}$ now refer to the regions left of the left shock wave and right of the right shock wave.
From conditions (\ref{EQU036})-(\ref{EQU035}) it can be shown that any analytic solution will have constant pressure $p^\ast$ and constant velocity $v^\ast$ between the two shock waves.
Let for $K\in\{L,R\}$ be
\begin{eqnarray}
P_{\mathrm{K}}(p,\rho_{\mathrm{K}}) &=& \sqrt{\frac{2(p-p_{\mathrm{K}})^2}{\rho_{\mathrm{K}}[\gamma(p+p_{\mathrm{K}})+(p-p_{\mathrm{K}})]}}
\end{eqnarray}
Then $p^\ast$ is implicit given by the root of the equation
\begin{eqnarray}
P_{\mathrm{L}}(p^\ast,\rho_{\mathrm{L}})+P_{\mathrm{R}}(p^\ast,\rho_{\mathrm{R}}) &=& u_{\mathrm{L}}-u_{\mathrm{R}}
\end{eqnarray}
and the velocity $v^\ast$ is given by
\begin{eqnarray}
v^\ast &=&\frac{v_{\mathrm{L}}+v_{\mathrm{R}}}{2}+\frac{P_{\mathrm{R}}(p^\ast,\rho_{\mathrm{R}})-P_{\mathrm{L}}(p^\ast,\rho_{\mathrm{L}})}{2}
\end{eqnarray}
Between the left and the right shock wave there is a contact wave separating the area between the two shock waves into \rm{III} and \rm{IV}.
% In contrary to the constant pressure $p^\ast=p_{\Rmnum{2}}=p_{\Rmnum{3}}$ and the constant velocity $v^\ast=v_{\Rmnum{2}}=v_{\Rmnum{3}}$ there is a discontinuity in $\rho$, which is only piecewise constant.
\begin{align}
&\rho_{\mathrm{III}}=\rho_{\mathrm{L}}\frac{\gamma(p^\ast+p_{\mathrm{L}})+(p^\ast-p_{\mathrm{L}})}{\gamma(p^\ast+p_{\mathrm{L}})-(p^\ast-p_{\mathrm{L}})}&
&\rho_{\mathrm{IV}} =\rho_{\mathrm{R}}\frac{\gamma(p^\ast+p_{\mathrm{R}})+(p^\ast-p_{\mathrm{R}})}{\gamma(p^\ast+p_{\mathrm{R}})-(p^\ast-p_{\mathrm{R}})}&\\
&v_{\mathrm{III}}=v^\ast&
&v_{\mathrm{IV}} =v^\ast&\\
&p_{\mathrm{III}}=p^\ast&
&p_{\mathrm{IV}} =p^\ast&
\end{align}
The propagation velocity of the left and the right shock waves are
\begin{eqnarray}
\sigma_{\mathrm{L}}&=&v_{\mathrm{L}}-\sqrt{\frac{\gamma(p^\ast+p_{\mathrm{L}})+(p^\ast-p_{\mathrm{L}})}{2\rho_{\mathrm{L}}}}\\
\sigma_{\mathrm{R}}&=&v_{\mathrm{R}}+\sqrt{\frac{\gamma(p^\ast+p_{\mathrm{R}})+(p^\ast-p_{\mathrm{R}})}{2\rho_{\mathrm{R}}}}
\end{eqnarray}
Taking advantage of the constant propagation speeds of the waves a similarity solution for $\frac{x}{t}$ can be given immediately
\begin{eqnarray}
\rho(t,x)=
\begin{cases}
\rho_{\mathrm{I}}=\rho_\mathrm{L} \qquad& x/t\le\sigma_{\mathrm{L}}\\
\rho_{\mathrm{III}} \qquad& \sigma_{\mathrm{L}}\le x/t\le v^\ast\\
\rho_{\mathrm{IV}} \qquad& v^\ast\le x/t\le \sigma_{\mathrm{R}}\\
\rho_{\mathrm{IV}}=\rho_\mathrm{R} \qquad& \sigma_{\mathrm{R}}\le x/t
\end{cases}
\end{eqnarray}
% Some final notes are appropriate at this point.
% The specific internal energy $U$ under the assumption of an ideal gas is only depended on the temperature $T$ of the ISM.
% In the simplest model without cooling and heating we can assume a constant temperature which, together with (\ref{EQU004}) results in a constant speed of sound $a_{\mathrm{L}}=a_{\mathrm{R}}=a$.
% However, the sound of speed $a$ is negligible small compared to the propagation speed $\sigma_{\mathrm{L}}$ and $\sigma_{\mathrm{R}}$ of the left and right shock wave.
% This pictures again that the corresponding eigenvalues $v\pm a$ to the shock waves are not their propagation speed.
For strong shock waves, the speed of sound $a$ is negligible compared to the propagation speeds $\sigma_{\mathrm{L}}$ and $\sigma_{\mathrm{R}}$
\begin{eqnarray}
a=\sqrt{\frac{\gamma p}{\rho}}=\sqrt{(\gamma-1)\gamma U}=\sqrt{\frac{\gamma k_\mathrm{B}T}{\mu_\mathrm{H}}}\ll\sigma_\mathrm{K}
\end{eqnarray}
and one can see clearly that for this nonlinear system of PDEs the eigenvalues $v\pm a$ are not the propagation speeds $\sigma_\mathrm{K}$ of the waves.
\subsubsection{Rarefaction solution}
The rarefaction solution is also based on $p^\ast$ and $v^\ast$, but is slightly more complicated.
Some differential schemes require the propagation velocities $\sigma$ for all kind of waves.
These are given in chapter \ref{WAF}.
For the full set of formulas we refer to \cite{TOR001}.
\subsection{Differential schemes}\label{SCHEMES}
In this chapter an overview of common differential schemes is given.
All differential schemes apply to a set of homogeneous PDEs.
If an additional non-homogeneous term is given, first the homogeneous system is forwarded and then the inhomogeneous part is treated as an ordinary differential equation.
PDEs with source terms depending on derivatives such as the Navier--Stokes equations will not be covered here, although there are possibilities to approximate quantities such as viscosity.
All differential schemes listed here follow the generic form
\begin{eqnarray}\label{EQU026}
\mathbf u_i^{(n+1)}=\mathbf u^{(n)}_i-\frac{\Delta t}{\Delta x}\left(\mathbf f_{i+\frac{1}{2}}-\mathbf f_{i-\frac{1}{2}}\right)
\end{eqnarray}
where the vectorial and scalar quantities build a staggered grid.
Fluxes with upper time-indices such as $\mathbf f^{(n)}_{i}=\mathbf f^{(n)}_{i}(\mathbf u^{(n)}_{i})$ are auxiliary quantities, while fluxes without these such as $\mathbf f^\mathrm{method}_{i+\frac{1}{2}}$ are final estimates to be used in (\ref{EQU026}).
\subsubsection{Lax--Friedrichs scheme}
\begin{eqnarray}
\mathbf f^{\text{LF}}_{i+\frac{1}{2}}=\frac{1}{2}\left(\mathbf f^{(n)}_{i}+\mathbf f^{(n)}_{i+1}\right)-\frac{\Delta x}{2\Delta t}\left(\mathbf u^{(n)}_{i+1}-\mathbf u^{(n)}_{i}\right)
\end{eqnarray}
In this particular scheme forwarding $\mathbf u$ turns out to be simple
\begin{eqnarray}
\mathbf u_i^{(n+1)}=\frac{1}{2}(\mathbf u_{i+1}^{(n)}+\mathbf u_{i-1}^{(n)})-\frac{\Delta t}{2\Delta x}\left(\mathbf f_{i+1}^{(n)}-\mathbf f_{i-1}^{(n)}\right)
\end{eqnarray}
For complex systems the numerical solutions are not very smooth and therefore the Lax--Friedrichs scheme is only used rarely.
\subsubsection{Lax--Wendroff--Richtmyer}
\begin{eqnarray}
\mathbf u_{i+\frac{1}{2}}^{(n+\frac{1}{2})}&=&\frac{1}{2}\left(\mathbf u^{(n)}_{i}+\mathbf u^{(n)}_{i+1}\right)-\frac{\Delta t}{2\Delta x}\left(\mathbf f^{(n)}_{i+1}-\mathbf f^{(n)}_{i}\right)\\
\mathbf f^{\text{LWR}}_{i+\frac{1}{2}}&=&\mathbf f_{i+\frac{1}{2}}^{(n+\frac{1}{2})}
\end{eqnarray}
This scheme is of second-order accuracy in space and time.
The drawback of its simplicity is a strong tendency to produce oscillations.
\subsubsection{First-order centred scheme (FORCE)}
\begin{eqnarray}
\mathbf f^{\text{FORCE}}_{i+\frac{1}{2}}&=&\frac{1}{2}\left(\mathbf f^{\text{LF}}_{i+\frac{1}{2}}+\mathbf f^{\text{LWR}}_{i+\frac{1}{2}}\right)
\end{eqnarray}
Here the fluxes are averaged over the fluxes of Lax--Friedrichs and Lax--Wendroff--Richtmyer.
This scheme lacks accuracy but is easy to implement and still useful for a large class of hyperbolic PDEs with smooth solutions.
\subsubsection{Lax--Wendroff--MacCormack}
This scheme is a Lax--Wendroff type scheme.
\begin{eqnarray}
\mathbf u_i^\ast &=& \mathbf u_j^{(n)}-\frac{\Delta t}{\Delta x}\left(\mathbf f_{i+1}^{(n)}-\mathbf f_{i}^{(n)}\right) \\
\mathbf f^{\text{LWMC}}_{i+\frac{1}{2}} &=& \frac{1}{2}\left(\mathbf f_{i+1}^{(n)}+\mathbf f_i^\ast\right)
\end{eqnarray}
In the linear case it becomes equivalent to the Lax--Wendroff--Richtmeyr scheme.
\subsubsection{Godunov scheme}
Let $\mathbf u_{i+\frac{1}{2}}^{(\text{RP})}(\frac{x}{t})$ be the solution\footnote{Since this is a self-similarity solution it is sufficient to use the argument $\frac{x}{t}$.} of the Riemann problem $\operatorname{RP}(\mathbf u_i^{(n)},\mathbf u_{i+1}^{(n)})$.
The approximated flux is its solution at $x=0$.
\begin{eqnarray}
\mathbf f^{\text{GOD}}_{i+\frac{1}{2}}&=&\mathbf f\left(\mathbf u_{i+\frac{1}{2}}^{(\text{RP})}(0)\right)
\end{eqnarray}
Here the calculation of the fluxes are considered as separated Riemann-problems.
This scheme is a monotone scheme and therefore only first-order accurate.
\subsubsection{Roe scheme}
This scheme was first introduced by Roe \cite{ROE001} in 1981 and is one of the most powerful schemes today.
It is a generalization of the Godunov scheme in the sense that the equations are again considered piecewise as Riemann problems.
The basic idea behind the Roe scheme is that the non-linear hyperbolic equations are first linearized.
As for the Riemann problem, an exact solution for the linear hyperbolic equations in general is known.
Better performance is achieved at the cost of greater complexity.
% The price for better performance is paid in form of higher complexity.
Roe's approach is to quasi-linearize $\partial_t\mathbf u+\vec\nabla f(\mathbf u)=0$ in the form $\partial_t\mathbf u+A\vec\nabla \mathbf u=0$ and give a good estimate for $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})$.
He postulated three conditions for his matrix $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})$:
\begin{enumerate}
\item $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})$ has only real eigenvalues with linear independent eigenvectors.
This is to ensure that the linearized equations are still hyperbolic.
\item $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})\cdot(\mathbf u_{\mathrm{R}}-\mathbf u_{\mathrm{L}})=f(\mathbf u_{\mathrm{R}})-f(\mathbf u_{\mathrm{L}})$.
This is to ensure conservation of quantities.
\item $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})\to \frac{\partial f(\mathbf u)}{\partial\mathbf u}$ as $\mathbf u_{\mathrm{R}}\to \mathbf u_{\mathrm{L}}$.
This is to ensure consistency.
\end{enumerate}
Remember that $\mathbf u=(\rho, \rho v, \epsilon)$ and that $p=(\gamma-1)(\epsilon-\frac{1}{2}\rho v^2)$.
It is furthermore convenient to use the total specific enthalpy $h=\frac{\epsilon+p}{\rho}$ and the adiabatic speed of sound $a=\sqrt{\frac{\gamma p}{\rho}}$.
We define the Roe-averages as
\begin{eqnarray}
\bar v &=& \frac{\sqrt{\rho_{\mathrm{L}}}v_{\mathrm{L}}+\sqrt{\rho_{\mathrm{R}}}v_{\mathrm{R}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}} = \frac{\sqrt{\rho_i^{(n)}}v_i^{(n)}+\sqrt{\rho_{i+1}^{(n)}}v_{i+1}^{(n)}}{\sqrt{\rho_{i}^{(n)}}+\sqrt{\rho_{i+1}^{(n)}}} \\
\bar h &=& \frac{\sqrt{\rho_{\mathrm{L}}}h_{\mathrm{L}}+\sqrt{\rho_{\mathrm{R}}}h_{\mathrm{R}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}} = \frac{\sqrt{\rho_i^{(n)}}h_i^{(n)}+\sqrt{\rho_{i+1}^{(n)}}h_{i+1}^{(n)}}{\sqrt{\rho_{i}^{(n)}}+\sqrt{\rho_{i+1}^{(n)}}} \\
\bar a &=& \sqrt{(\gamma-1)\left(\bar h-\frac{\bar v^2}{2}\right)} \\
\bar \rho &=& \sqrt{\rho_{\mathrm{L}}\rho_{\mathrm{R}}} = \sqrt{\rho_{i}^{(n)}\rho_{i+1}^{(n)}}
\end{eqnarray}
For the Euler equations as given in (\ref{EQU001})-(\ref{EQU003}) the matrix $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})$ is uniquely determined by Roe's postulated requirements:
\begin{eqnarray}
A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}}) &=&
\begin{pmatrix}
0 & 1 & 0 \\
\frac{\gamma-3}{2}\bar v^2 & (3-\gamma)\bar v & \gamma-1 \\
\frac{\gamma-3}{2}\bar v^3-\bar v\bar h & \bar h-(\gamma-1)\bar v^2 & \gamma\bar v &
\end{pmatrix}
\end{eqnarray}
The eigenvalues of $A(\mathbf u_{\mathrm{L}},\mathbf u_{\mathrm{R}})$ are
\begin{eqnarray}
\mathbf\lambda &=&
\begin{pmatrix}
\bar v-\bar a \\
\bar v \\
\bar v+\bar a
\end{pmatrix}
\end{eqnarray}
and the corresponding matrix of eigenvectors is
\begin{eqnarray}
R &=&
\begin{pmatrix}
1 & 1 & 1 \\
\bar v-\bar a & \bar v & \bar v+\bar a \\
\bar h-\bar v\bar a & \frac{\bar v^2}{2} & \bar h+\bar v\bar a
\end{pmatrix}
\end{eqnarray}
We define a difference vector
\begin{eqnarray}
\mathbf d &=& \frac{1}{\bar a^2}
\begin{pmatrix}
\frac{1}{2}(\Delta p-\bar a\bar\rho\Delta v) \\
\bar a^2\Delta\rho-\Delta p \\
\frac{1}{2}(\Delta p+\bar a\bar\rho\Delta v)
\end{pmatrix}
\end{eqnarray}
with $\Delta \rho=\rho_{\mathrm{R}}-\rho_{\mathrm{L}}$, $\Delta v=v_{\mathrm{R}}-v_{\mathrm{L}}$ and $\Delta p=p_{\mathrm{R}}-p_{\mathrm{L}}$ as usual and finally multiply each component with the absolute value of the corresponding eigenvalue.
\begin{eqnarray}
\mathbf d^\ast &=&
\begin{pmatrix}
|\lambda_1| d_1 \\
|\lambda_2| d_2 \\
|\lambda_3| d_3 \\
\end{pmatrix}
\end{eqnarray}
The corrected Roe-flux for the Roe scheme reads
\begin{eqnarray}
\mathbf f^{\mathrm{ROE}}_{i+\frac{1}{2}} &=& \frac{1}{2}\left(\mathbf f^{(n)}_{i+1}+\mathbf f^{(n)}_{i}\right)-\frac{1}{2}\,R\cdot\mathbf d^\ast
\end{eqnarray}
\subsubsection{Steger--Warming scheme}
A scheme similar to the Roe scheme but easier to implement is the Steger--Warming scheme.
This scheme splits the flux into a positive and a negative part by testing the signatures of the eigenvalues of $\frac{\partial f(\mathbf u)}{\partial \mathbf u}$.
Let
\begin{eqnarray}
R(\mathbf u) &=&
\begin{pmatrix}
1 & 1 & 1 \\
v-a & v & v+a \\
h-v a & \frac{v^2}{2} & h+v a
\end{pmatrix}
\end{eqnarray}
with total specific enthalpy $h=\frac{\epsilon+p}{\rho}$ and the vector of eigenvalues $\mathbf\lambda=(v-a,v,v+a)^{T}$.
We now distinguish between the positive and the negative eigenvalues by defining
\begin{eqnarray}
\lambda^\pm = \frac{1}{2}(\lambda\pm |\lambda|)
\end{eqnarray}
We further define two vectors $\mathbf d^\pm$ by
\begin{eqnarray}
\mathbf d^\pm(\mathbf u) &=& \frac{\rho}{2\gamma}
\begin{pmatrix}
\lambda_1^\pm \\
2(\gamma-1)\lambda_2^\pm \\
\lambda_3^\pm \\
\end{pmatrix}
\end{eqnarray}
and the separated fluxes by
\begin{eqnarray}
\mathbf f^\pm &=& R\cdot\mathbf d^\pm
\end{eqnarray}
The Steger--Warming flux reads
\begin{eqnarray}
\mathbf f^{\mathrm{SW}}_{i+\frac{1}{2}} &=& \mathbf {f^+}^{(n)}_{i}+\mathbf {f^-}^{(n)}_{i+1}
\end{eqnarray}
\subsubsection{Van-Leer scheme}\label{SCHEMESVL}
As in the Steger--Warming scheme the flux vector is split in an upper and a lower part
\begin{eqnarray}
\mathbf f^\pm &=& \pm\frac{\rho a}{4}\left(1\pm\frac{v}{a}\right)^2
\begin{pmatrix}
1 \\
\pm\frac{2a}{\gamma}(1\pm\frac{\gamma-1}{2}\frac{v}{a}) \\
\frac{2a^2}{\gamma^2-1}(1\pm\frac{\gamma-1}{2}\frac{v}{a})^2
\end{pmatrix}
\end{eqnarray}
and the Van Leer flux is
\begin{eqnarray}
\mathbf f^{\mathrm{VL}}_{i+\frac{1}{2}} &=& \mathbf {f^+}^{(n)}_{i}+\mathbf {f^-}^{(n)}_{i+1}
\end{eqnarray}
\subsubsection{Weighted average flux (WAF) scheme}\label{WAF}
Let $\mathbf u_{i+\frac{1}{2}}^{(\text{RP})}(\frac{x}{t})$ be the solution of the Riemann problem $\operatorname{RP}(\mathbf u_i^{(n)},\mathbf u_{i+1}^{(n)})$ as explained in chapter \ref{RIEMANNPROBLEM}.
Remember that there are three waves partitioning the solution domain into four regions.
In addition, the left and/or right wave might be a rarefaction wave, which further partitions the solution.
This gives from left to right up to six regions for each of which the solution $\mathbf u_{i+\frac{1}{2}}^{(\text{RP})}(\frac{x}{t})$ can be calculated with arbitrary precision.
We will denote the boundaries between these regions by $\omega={1,2,3,4,5}$ and note that in the case of a left (right) shock wave the boundaries 1 and 2 (4 and 5) overlap.
Let $\sigma^{(\omega)}$ be the propagation velocities of $\omega$ (see chapter \ref{RIEMANNPROBLEM})
\begin{eqnarray}
\sigma^{(1)} &=&
\begin{cases}
v_\mathrm{L}-a_\mathrm{L}\sqrt{\frac{p^\ast}{p_\mathrm{L}}\frac{\gamma+1}{2\gamma}+\frac{\gamma-1}{2\gamma}} & \text{shock wave} \\
v_\mathrm{L}-a_\mathrm{L} & \text{rarefaction wave}
\end{cases} \\
\sigma^{(2)} &=&
\begin{cases}
\sigma^{(1)} & \,\,\,\,\qquad\text{shock wave} \\
v_\mathrm{L}-a_\mathrm{L}\left(\frac{p^\ast}{p_\mathrm{L}}\right)^{\frac{\gamma-1}{2\gamma}} &\,\,\,\,\qquad\text{rarefaction wave}
\end{cases} \\
\sigma^{(3)} &=& v^\ast\\
\sigma^{(4)} &=&
\begin{cases}
\sigma^{(5)} & \,\,\,\,\qquad\text{shock wave} \\
v_\mathrm{R}+a_\mathrm{R}\left(\frac{p^\ast}{p_\mathrm{R}}\right)^{\frac{\gamma-1}{2\gamma}} & \,\,\,\,\qquad\text{rarefaction wave}
\end{cases} \\
\sigma^{(5)} &=&
\begin{cases}
v_\mathrm{R}+a_\mathrm{R}\sqrt{\frac{p^\ast}{p_\mathrm{R}}\frac{\gamma+1}{2\gamma}+\frac{\gamma-1}{2\gamma}} & \text{shock wave} \\
v_\mathrm{R}+a_\mathrm{R} & \text{rarefaction wave}
\end{cases}
\end{eqnarray}
with $p^\ast$ and $v^\ast$ from the solution given in chapter \ref{RIEMANNPROBLEM} and $a=\sqrt{\frac{\gamma p}{\rho}}$ the adiabatic speed of sound as usually.
The Courant number $c^{(\omega)}$ for each of these velocities is defined by
\begin{eqnarray}\label{EQU027}
c^{(\omega)} &=& \frac{\Delta t}{\Delta x}\sigma^{(\omega)}
\end{eqnarray}
For a quantity $q$, let $\Delta^{(\omega)}q$ be the step that $q$ takes at $\omega$, thus e.g.\ $\Delta^{(\omega)}\mathbf f=\mathbf f^{(\omega)}_\mathrm{R}-\mathbf f^{(\omega)}_\mathrm{L}$.
The weighted average flux then reads as follows
\begin{eqnarray}\label{EQU031}
\mathbf f^{\mathrm{(WAF)}}_{i+\frac{1}{2}} &=&
\frac{1}{2}\left(\mathbf f_i^{(n)}+\mathbf f_{i+1}^{(n)}\right)-\frac{1}{2}\sum_\omega c^{(\omega)}_{i+\frac{1}{2}}\Delta^{(\omega)}\mathbf f^{\mathrm{(RP)}}_{i+\frac{1}{2}}
\end{eqnarray}
A small variation is possible by averaging the non-conservative set $\mathbf w=(\rho, v, p)$ instead of the fluxes
\begin{eqnarray}\label{EQU032}
\mathbf w^{\mathrm{(WAF)}}_{i+\frac{1}{2}} &=&
\frac{1}{2}\left(\mathbf w_i^{(n)}+\mathbf w_{i+1}^{(n)}\right)-\frac{1}{2}\sum_\omega c^{(\omega)}_{i+\frac{1}{2}}\Delta^{(\omega)}\mathbf w^{\mathrm{(RP)}}_{i+\frac{1}{2}} \\
\mathbf f^{\mathrm{(WAF)}}_{i+\frac{1}{2}} &=& f(\mathbf w^{\mathrm{(TVD)}}_{i+\frac{1}{2}})
\end{eqnarray}
\subsubsection{Total variation diminishing (TVD)}\label{SCHEMESTVD}
Non-monotonic schemes will sooner or later cause oscillations in the solutions.
TVD-schemes attack this problem by ensuring that the total variation of the numerical solution within a cell is not increasing.
This can be accomplished with so-called flux limiter functions.
In the following a TVD method for the WAF scheme is described together with some possible flux limiters.
Let again $\Delta^{(\omega)}q$ be the step that $q$ takes at $\omega$, thus $q^{(\omega)}_{\mathrm{R}}-q^{(\omega)}_{\mathrm{L}}$.
For $q=\rho^{\mathrm{(RP)}}$ we define the flow parameter $r^{(\omega)}$ by
\begin{eqnarray}\label{EQU028}
r^{(\omega)}_{i+\frac{1}{2}}=r^{(\omega)}_{i+\frac{1}{2},\rho} &=&
\begin{cases}
\frac{\Delta^{(\omega)}q_{i+\frac{3}{2}}}{\Delta^{(\omega)}q_{i+\frac{1}{2}}}=
\frac{\Delta^{(\omega)}\rho_{i+\frac{3}{2}}^{\mathrm{(RP)}}}{\Delta^{(\omega)}\rho_{i+\frac{1}{2}}^{\mathrm{(RP)}}}
& \text{for }c^{(\omega)}<0 \\
\frac{\Delta^{(\omega)}q_{i-\frac{1}{2}}}{\Delta^{(\omega)}q_{i+\frac{1}{2}}}=
\frac{\Delta^{(\omega)}\rho_{i-\frac{1}{2}}^{\mathrm{(RP)}}}{\Delta^{(\omega)}\rho_{i+\frac{1}{2}}^{\mathrm{(RP)}}}
& \text{for }c^{(\omega)}>0
\end{cases}
\end{eqnarray}
Alternatively, one can also choose $q=\epsilon^{\mathrm{(RP)}}$ to get $r^{(\omega)}_{i+\frac{1}{2},\epsilon}$.
The flow parameter $r^{(\omega)}$, together with the Courant number $c^{(\omega)}$ introduced in chapter \ref{WAF} serves as an input for WAV flux limiter functions $\phi(r,c)$.
This suffices for the calculation of the TVD flux
\begin{eqnarray}\label{EQU029}
\mathbf f^{\mathrm{(TVD)}}_{i+\frac{1}{2}} &=&
\frac{1}{2}\left(\mathbf f_i^{(n)}+\mathbf f_{i+1}^{(n)}\right)-\frac{1}{2}\sum_\omega\operatorname{sign}(c^{(\omega)}_{i+\frac{1}{2}})\phi(r^{(\omega)}_{i+\frac{1}{2}},c^{(\omega)}_{i+\frac{1}{2}})\Delta^{(\omega)}\mathbf f^{\mathrm{(RP)}}_{i+\frac{1}{2}}
\end{eqnarray}
As for the WAF-scheme, a small variation is possible by averaging the non-conservative set $\mathbf w=(\rho, v, p)$ instead of the fluxes
\begin{eqnarray}\label{EQU030}
\mathbf w^{\mathrm{(TVD)}}_{i+\frac{1}{2}} &=&
\frac{1}{2}\left(\mathbf w_i^{(n)}+\mathbf w_{i+1}^{(n)}\right)-\frac{1}{2}\sum_\omega\operatorname{sign}(c^{(\omega)}_{i+\frac{1}{2}})\phi(r^{(\omega)}_{i+\frac{1}{2}},c^{(\omega)}_{i+\frac{1}{2}})\Delta^{(\omega)}\mathbf w^{\mathrm{(RP)}}_{i+\frac{1}{2}} \\
\mathbf f^{\mathrm{(TVD)}}_{i+\frac{1}{2}} &=& f(\mathbf w^{\mathrm{(TVD)}}_{i+\frac{1}{2}})
\end{eqnarray}
Note that by replacing the term $\operatorname{sign}(c^{(\omega)}_{i+\frac{1}{2}})\phi(r^{(\omega)}_{i+\frac{1}{2}},c^{(\omega)}_{i+\frac{1}{2}})$ in the equations (\ref{EQU029}) or (\ref{EQU030}) by $c^{(\omega)}_{i+\frac{1}{2}}$ one obtains the simpler WAF-scheme from above.
The WAV flux limiter functions $\phi(r,c)$ are correlated with conventional flux limiters $\psi(r)$ via
\begin{eqnarray}
\phi(r,c) &=& 1-(1-|c|)\psi(r)
\end{eqnarray}
% There are as many as 15 flux limiters $\psi(r)$ listed in Wikipedia and probably there are even more.
There is a minimum of 15 feasible flux limiters $\psi(r)$ with many listed in \cite{TOR001}.
Two commonly used examples are $\mathrm{SUPERBEE}$ (sa) and $\mathrm{van\,Leer}$ (vl).
\begin{eqnarray}
\psi_{\mathrm{sa}}(r) &=&\max\left[0,\min\left(2r,1\right),\min\left(r,2\right)\right] \\
\psi_{\mathrm{vl}}(r) &=&\frac{r+\left|r\right|}{1+\left|r\right|}
\end{eqnarray}
% (Toro, p456, p496)
\subsubsection{Other schemes}
There is a long list of additional schemes including:
\begin{itemize}
\item Monotone upstream-centred schemes for conservation laws (MUSCL)
\item Osher scheme
\item Kurganov and Tadmor central scheme
\item Liou-Steffen advection upstream splitting method (AUSM) scheme.
\end{itemize}
and probably many others.
From this list, the MUSCL schemes are probably the most relevant ones.
% (Toro, p504)
There is a large number of variations of MUSCL schemes.
The most common is the MUSCL-Hancock Method (MHM), which is exhaustively explained in \cite{TOR001}.
It is a rather complex three step second-order accurate scheme and promises excellent performance for further investigations.
Table \ref{TAB001} gives a summary of the schemes introduced here.
\begin{table}
\begin{center}
\begin{tabular}{lll}
Acronym & Scheme & Variation \\
\hline
LF & Lax--Friedrichs & \\
LW-R & Lax--Wendroff & Richtmeyr \\
LW-MC & Lax--Wendroff & MacCormack \\
FORCE & First-Order Centred Scheme & \\
GOD & Godunov & \\
ROE & Roe & \\
SW & Steger--Warming & \\
VL & van Leer & \\
WAF-F & Weighted average flux & flux $\mathbf f$ averaged \\
WAF-W & Weighted average flux & $\mathbf w$ averaged \\
TVD-F-SA & Total variation diminishing & flux $\mathbf f$ averaged, SUPERBEE limiter \\
TVD-F-VL & Total variation diminishing & flux $\mathbf f$ averaged, van Leer limiter \\
TVD-W-SA & Total variation diminishing & $\mathbf w$ averaged, SUPERBEE limiter \\
TVD-W-VL & Total variation diminishing & $\mathbf w$ averaged, van Leer limiter
\end{tabular}
\caption{List of schemes under investigation. The TVD-schemes are based on the WAF approach.}\label{TAB001}
\end{center}
\end{table}
\subsection{Self-similar solution}\label{SEDOVSOLUTION}
The self-similar solution of an intense explosion was comprehensively discussed in 1959 by Sedov in the first edition of \cite{SED002}.
His work is based on the textbook of Landau \& Lifschitz about Fluid Mechanics \cite{LAN001} and develops solutions for a class of exponentially decreasing initial density distributions.
Nevertheless, a solution was already given in 1950 by Taylor \cite{TAY001}\cite{TAY002} in the context of nuclear bomb blasts, but has not been adapted for supernovae.
The derivation of the exact self-similarity solution is cumbersome and needs a lot of work.
An excellent summary of the results is given by \cite{KAM001}, but even this covers about 16 pages.
I will therefore stick to two special cases, the one with constant initial density distribution and the one with exponentially decreasing distribution with exponent $\omega=-1.5$ and will give solutions for them.
In both cases I assume an ideal monoatomic gas with $\gamma=\frac{5}{3}$.
Fortunately, with a given initial density distribution
\begin{eqnarray}
\rho_0 &=& A\,r^{-\omega}
\end{eqnarray}
together with the energy of the blast $E^\ast$ the solution is already uniquely determined.
We will choose the constant $A$ consistent with the mass $M^\ast$ of the supernova and the density $\rho_\mathrm{ISM}$ of the interstellar medium.
The exponent $\omega$ cannot be chosen arbitrarily, and actually many values will cause problems such as infinite mass, vacuum solutions or other singularities.
We note that $\omega<2$ will work for $\gamma=\frac{5}{3}$ and refer for details to \cite{KAM001}.
The solution at some time $t$ will usually have the parameterized form $r(\lambda),\,v(\lambda),\,\rho(\lambda),\,p(\lambda)$, where $\lambda=0$ corresponds to the origin and $\lambda=1$ corresponds to the shell of the shock wave.
Neither arguments $\lambda>1$ nor $\lambda<0$ make sense.
Since our major objective of deriving an exact self-similarity solution is to get useful initial values for the simulation of the Sedov-Taylor phase we have to depart from the common approach.
We will initially give the position of the shell $r_\mathrm{shell}$ and evaluate the time $t_0$ at which the shell is at $r_\mathrm{shell}$.
The other physical quantities will follow from this.
However, this is not a trivial task and requires the use of the Rankine-Hugoniot equations as well as the numerical evaluation of two integrals, where, on top of everything, at least one has an infinite value at one limit of integration.
\subsubsection{Constant initial density}
The solution for constant initial density $\rho_0\equiv A$ and specified blast energy $E^\ast$ is given by
\begin{eqnarray}
% \frac{r(V)}{r_\mathrm{shell}} &=& 2^{-\frac{16}{65}}\cdot 3^{\frac{16}{65}}\cdot 5^{-\frac{20}{39}}\cdot V^{-\frac{2}{5}}\cdot \left(V-\frac{6}{25}\right)^{\frac{2}{13}}\cdot \left(\frac{1}{2}-V\right)^{-\frac{82}{195}} \\
% \approx\\
% &\approx& 0.484\cdot V^{-\frac{2}{5}}\cdot \left(V-0.24\right)^{\frac{2}{13}}\cdot \left(0.5-V\right)^{-\frac{82}{195}}\\
% \frac{v(V)}{v_\mathrm{behind}} &=& 2^{\frac{49}{65}}\cdot 3^{-\frac{49}{65}}\cdot 5^{\frac{19}{39}}\cdot V^{\frac{3}{5}}\cdot \left(V-\frac{6}{25}\right)^{\frac{2}{13}}\cdot \left(\frac{1}{2}-V\right)^{-\frac{82}{195}} \\
% \frac{\rho(V)}{\rho_\mathrm{behind}} &=& 2^{-\frac{69}{13}}\cdot 3^{-\frac{9}{13}}\cdot 5^{\frac{22}{13}}\cdot\left(V-\frac{6}{25}\right)^{\frac{9}{13}}\cdot \left(\frac{1}{2}-V\right)^{-\frac{82}{13}}\cdot \left(\frac{2}{5}-V\right)^{-6} \\
\frac{r(\lambda)}{r_\mathrm{shell}} &=& \lambda^{\frac{2}{13}}\cdot\left(\frac{13-3\lambda}{10}\right)^{-\frac{82}{195}}\cdot\left(\frac{\lambda+4}{5}\right)^{-\frac{2}{5}} \\
\frac{v(\lambda)}{v_\mathrm{behind}} &=& \lambda^{\frac{2}{13}}\cdot\left(\frac{13-3\lambda}{10}\right)^{-\frac{82}{195}}\cdot\left(\frac{\lambda+4}{5}\right)^{\frac{3}{5}} \\
\frac{\rho(\lambda)}{\rho_\mathrm{behind}} &=& \lambda^{\frac{9}{13}}\cdot\left(\frac{13-3\lambda}{10}\right)^{\frac{82}{13}}\cdot\left(\frac{8-3\lambda}{5}\right)^{-6} \\
\frac{p(\lambda)}{p_\mathrm{behind}} &=& \left(\frac{13-3\lambda}{10}\right)^{\frac{82}{15}}\cdot\left(\frac{\lambda+4}{5}\right)^{\frac{6}{5}}\cdot\left(\frac{8-3\lambda}{5}\right)^{-5}
\end{eqnarray}
with
\begin{align}
& r_\mathrm{shell}=\left(\frac{E^\ast}{\alpha A}\right)^\frac{1}{5}\cdot t^{\frac{2}{5}}, &
& v_\mathrm{behind}=\frac{3}{10}\left(\frac{E^\ast}{\alpha A}\right)^\frac{1}{5}\cdot t^{-\frac{3}{5}}, &
& \rho_\mathrm{behind}=4 A, &
& p_\mathrm{behind}=\frac{3}{25}A\left(\frac{E^\ast}{\alpha A}\right)^\frac{2}{5}\cdot t^{-\frac{6}{5}} &
\end{align}
and
\begin{eqnarray}
\alpha=0.4935901495825678645\dots
\end{eqnarray}
Calculation of $\alpha$ is tricky and can only be done numerically.
Detailed instructions are given e.g.\ in \cite{KAM001}.
% This solution fulfills the Euler equations (\ref{EQU001})-(\ref{EQU003})\footnote{If the reader has the same enthusiasms about verifying the solution as the author, he will have to use $\frac{dr}{dt}=\frac{2}{5}\frac{r}{t}, \frac{dv}{dt}=-\frac{3}{5}\frac{v}{t}, \frac{d\rho}{dt}=0$.
This solution fulfills the Euler equations (\ref{EQU001})-(\ref{EQU003})\footnote{If the reader is as keen on verifying the solution as the author, he will have to use $\frac{dr}{dt}=\frac{2}{5}\frac{r}{t}, \frac{dv}{dt}=-\frac{3}{5}\frac{v}{t}, \frac{d\rho}{dt}=0$.
With this one can rewrite the partial time derivatives as $\frac{\partial\rho}{\partial t}=-\frac{\partial\rho}{\partial r}\frac{2}{5}\frac{r}{t}$ and $\frac{\partial v}{\partial t}=-\frac{3}{5}\frac{v}{t}-\frac{\partial v}{\partial r}\frac{2}{5}\frac{r}{t}$.
A good symbolical calculation software is extremely helpful.}.
Figure \ref{FIG010} illustrates the full self-similarity solution for a constant initial density distribution.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm]{SN.fig010.eps}
\includegraphics[width=8cm]{SN.fig011.eps}\\
\includegraphics[width=8cm]{SN.fig012.eps}
\includegraphics[width=8cm]{SN.fig013.eps}\\
\caption{$\lambda$-normalized self-similarity solution with density $\rho^\ast$, expansion velocity behind the shell $v^\ast$ and pressure $p^\ast$ for a constant initial density distribution $\rho^\ast_\mathrm{ISM}=0.0153$ at $t^\ast_0=0.0295$. The shell is currently located at $r^\ast=3.604 [\mathrm{pc}]$.}\label{FIG010}
\end{center}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm]{SN.fig014.eps}
\includegraphics[width=8cm]{SN.fig015.eps}\\
\includegraphics[width=8cm]{SN.fig016.eps}
\includegraphics[width=8cm]{SN.fig017.eps}\\
\caption{$\lambda$-normalized self-similarity solution with density $\rho^\ast$, expansion velocity behind the shell $v^\ast$ and pressure $p^\ast$ for an exponentially decreasing initial density distribution ($\omega=-1.5$) at $t^\ast_0=0.0245$. The shell is currently located at $r^\ast=2.860 [\mathrm{pc}]$.}\label{FIG014}
\end{center}
\end{figure}
In some of the simulations in chapter \ref{SEC005} we will use this solution as initial condition at $t_0$.
We assume that the total mass within $r_\mathrm{shell}$ at $t_0$ is $M^\ast=3$ and that the mass is equally distributed with density $\rho_\mathrm{ISM}=1\,\mathrm{cm}^{-3}$.
So there is no initial discontinuity in the density $\rho$.
The blast energy should be $E=10^{44}\,\mathrm{J}$.
This gives initial values at $t_0$ ($t_4=10^4\,\mathrm{years}$):
\begin{eqnarray}
t^\ast_0 &=& 0.0295497\,[t_4]\qquad(\cong 295\,\text{years})\\
r^\ast_\mathrm{shell} &=& 3.60396\,[\mathrm{pc}]\qquad(\cong 11.75\,\text{light-years})\\
v^\ast_\mathrm{behind} &=& 36.5888\,[\mathrm{pc}/t_4]\qquad(\cong 3578\,\mathrm{km}/\mathrm{s}) \\
\rho^\ast_\mathrm{behind} &=& 0.0612224\,[M_\odot/\mathrm{pc}^3]\qquad(\cong 4.5\,\mathrm{kg}/\text{volume of earth}) \\
p^\ast_\mathrm{behind} &=& 27.3203\,[M_\odot/(\mathrm{t_4^2\,pc})]\qquad(\cong 0.018\,\mu\mathrm{Pa})
\end{eqnarray}
\subsubsection{Decreasing initial density}
Alternatively we now consider an exponentially decreasing initial density profile
\begin{eqnarray}
\rho_0(r) &=&
\begin{cases}
A\,r^{-1.5} & \text{if } r \le r_\mathrm{shell} \\
\rho_\mathrm{ISM} & \text{if } r > r_\mathrm{shell} \\
\end{cases}
\end{eqnarray}
where
\begin{eqnarray}
A = \sqrt{\frac{3M\rho_\mathrm{ISM}}{8\pi}}
\qquad\text{and}\qquad
r_\mathrm{shell} =\sqrt[3]{\frac{3M}{8\pi\rho_\mathrm{ISM}}}
\end{eqnarray}
follow from the condition of a continuous density function and a total mass $M$ within $r_\mathrm{shell}$.
In this case, the self-similarity solution for specified blast energy $E^\ast$ is given by
\begin{eqnarray}
\frac{r(\lambda)}{r_\mathrm{shell}} &=& \lambda^{\frac{4}{11}}\cdot\left(\frac{11-6\lambda}{5}\right)^{-\frac{139}{231}}\cdot\left(\frac{\lambda+4}{5}\right)^{-\frac{4}{7}} \\
\frac{v(\lambda)}{v_\mathrm{behind}} &=& \lambda^{\frac{4}{11}}\cdot\left(\frac{11-6\lambda}{5}\right)^{-\frac{139}{231}}\cdot\left(\frac{\lambda+4}{5}\right)^{\frac{3}{7}} \\
\frac{\rho(\lambda)}{\rho_\mathrm{behind}} &=& \lambda^{\frac{3}{11}}\cdot\left(\frac{11-6\lambda}{5}\right)^{-\frac{417}{77}}\cdot\left(\frac{\lambda+4}{5}\right)^{\frac{6}{7}}\cdot\left(\frac{8-3\lambda}{5}\right)^{4} \\
\frac{p(\lambda)}{p_\mathrm{behind}} &=& \left(\frac{11-6\lambda}{5}\right)^{-\frac{139}{21}}\cdot\left(\frac{\lambda+4}{5}\right)^{\frac{12}{7}}\cdot\left(\frac{8-3\lambda}{5}\right)^{5}
\end{eqnarray}
with
\begin{align}
& r_\mathrm{shell}=\left(\frac{E^\ast}{\alpha A}\right)^\frac{2}{7}\cdot t^{\frac{4}{7}}, &
& v_\mathrm{behind}=\frac{3}{7}\left(\frac{E^\ast}{\alpha A}\right)^\frac{2}{7}\cdot t^{-\frac{3}{7}}, &
& \rho_\mathrm{behind}=4 A, &
& p_\mathrm{behind}=\frac{3}{49}A\left(\frac{E^\ast}{\alpha A}\right)^\frac{4}{7}\cdot t^{-\frac{6}{7}} &
\end{align}
and
\begin{eqnarray}
\alpha=1.07901668981107738115\dots
\end{eqnarray}
Using again $M^\ast=3$, $E=10^{44}\,\mathrm{J}$ and $\rho_\mathrm{ISM}=1\,\mathrm{cm}^{-3}$ this gives
\begin{eqnarray}
t^\ast_0 &=& 0.0245173\,[t_4]\qquad(\cong 245\,\text{years})\\
r^\ast_\mathrm{shell} &=& 2.86012\,[\mathrm{pc}]\qquad(\cong 9.33\,\text{light-years})\\
v^\ast_\mathrm{behind} &=& 49.9959\,[\mathrm{pc}/t_4]\qquad(\cong 4889\,\mathrm{km}/\mathrm{s}) \\
\rho^\ast_\mathrm{behind} &=& 0.0612224\,[M_\odot/\mathrm{pc}^3]\qquad(\cong 4.5\,\mathrm{kg}/\text{volume of earth}) \\
p^\ast_\mathrm{behind} &=& 51.0103\,[M_\odot/(\mathrm{t_4^2\,pc})]\qquad(\cong 0.033\,\mu\mathrm{Pa})
\end{eqnarray}
This solution is illustrated in figure \ref{FIG014}.
\newpage
.
\newpage
\section{A test case for numerical schemes}\label{SEC004}
\begin{abstract}
{\it
This chapter deals with a simple test case to study and verify the differential schemes introduced in chapter \ref{SCHEMES}.
This test case will also have the function of a sieve to eliminate methods that are not usable for the Euler equations and hence for the supernova model.
}
\end{abstract}
\subsection{Challenging the schemes}
It was a deliberate decision to include this chapter in this thesis since more fundamental astrophysical questions such as instabilities (e.g.\ Kelvin-Helmholtz instability), composition of the ISM, metallicity, agglomeration (e.g.\ collision of waves) and so on are only touched upon very perfunctorily.
The decisive factor was the complexity of the supernova model.
Without several years of experience in this field, it is hard to judge and interpret the results.
There are many stumbling blocks and things can and will go wrong during the first few attempts.
Many failures are not even apparent and will cause subsequent errors.
Before starting with complex supernovae simulations we will therefore insert some general investigations on differential schemes.
The model under investigation is the simple Riemann problem
\begin{align}
\rho_{\mathrm{L}}=1 &\qquad \rho_{\mathrm{R}}=0.125 &\qquad \text{initial density}\\
v_{\mathrm{L}}=0 &\qquad v_{\mathrm{R}}=0 &\qquad \text{initial velocity}\\
p_{\mathrm{L}}=1 &\qquad p_{\mathrm{R}}=0.1 &\qquad \text{initial pressure}\\
\end{align}
The exact solution is known and the quality of the results can therefore be quantified easily.
In this model $N=1001$ equidistant spatial points in the interval $[-1,1]$ are used with the initial discontinuity at the origin, thus $\Delta x=0.02$.
We will focus our interest to the numerical solution at $t=0.5$.
For all 14 schemes in table \ref{TAB001} the Courant-Friedrichs-Lewy number (see chapter \ref{SECCFL}) is set to $\operatorname{CFL}=0.99$.
In figure \ref{FIG018} and \ref{FIG019}, solutions for density $\rho$, velocity $v$, pressure $p$ and specific internal energy $\epsilon$ and some zooms for critical ranges are plotted.
It is interesting to observe that, as an effect of conservation, the oscillations for $\rho$, $p$ and $\epsilon$ follow the same pattern.
Basically we can observe two trends:
\begin{itemize}
\item smooth solutions with poor accuracy
\item oscillating solutions with high accuracy
\end{itemize}
To quantify the quality of the solution, we introduce an error measure $\varepsilon_\alpha$ and apply it to the different density solutions.
Let $\Delta\rho_i$ be the difference between the calculated and the exact solution at position $i\in\{0,1,\dots N-1\}$.
We define
\begin{eqnarray}
\varepsilon^{(\alpha)}_\rho &=&
\begin{cases}
\sqrt[\alpha]{\frac{\sum |\Delta\rho_i|^\alpha}{N}} & \qquad \alpha<\infty \\
\max(\Delta\rho_i) & \qquad \alpha=\infty
\end{cases}
\end{eqnarray}
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig018.Test1.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig019.Test1.v.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig020.Test1.p.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig021.Test1.e.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig022.Test1.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig023.Test1.v.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig024.Test1.p.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig025.Test1.e.eps}
\end{tabular}
\caption{Solutions for the test case of a simple Riemann problem: low-order schemes are smooth, but inaccurate, high-order scheme cause oscillations.
Cases of qualitatively similar solutions are indicated by CI, CII and CIII.}\label{FIG018}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig026.Test1.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig027.Test1.v.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig030.Test1.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig031.Test1.v.eps} \\
\end{tabular}
\caption{More magnifications of solutions for the test case of a simple Riemann problem.}\label{FIG019}
\end{center}
\end{figure}
\begin{table}[H]
\begin{center}
\begin{tabular}{r|ccc|l}
Scheme & $N\epsilon^{(1)}_\rho$ & $N\epsilon^{(2)}_\rho$ & $N\epsilon^{(\infty)}_\rho$ & problem \\
\hline
LF & 9.339 & 20.640 & 129.476 & not smooth, curled \\
LW-R & 2.223 & 9.713 & 143.920 & heavy oscillations \\
LW-MC & 2.273 & 9.883 & 140.222 & heavy oscillations \\
FORCE & 6.114 & 16.612 & 125.237 & low accuracy \\
GOD & 4.205 & 13.825 & 121.416 & low accuracy \\
ROE & 4.219 & 13.785 & 122.363 & oscillations \\
SW & 5.740 & 15.811 & 131.276 & oscillations \\
VL & 4.623 & 14.167 & 125.896 & others \\
WAF-F & 2.134 & 9.728 & 140.633 & oscillations \\
WAF-W & 2.120 & 9.700 & 146.688 & oscillations \\
TVD-F-SA & 0.511 & 5.304 & 129.605 & complex \\
TVD-F-VL & 0.911 & 7.082 & 129.716 & complex \\
TVD-W-SA & 0.591 & 4.530 & 110.865 & complex \\
TVD-W-VL & 0.990 & 6.521 & 118.375 & complex
\end{tabular}
\end{center}
\caption{Quality of different schemes for the test case of a simple Riemann problem. The factor $N$ is used to get convenient numbers.}\label{TAB002}
\end{table}
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig032.Test1.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig033.Test1.v.eps} \\
\end{tabular}
\caption{A very close magnification of the test case solution for TVD-F-SA. The oscillations are reduced by using the complex approach introduced in chapter \ref{SCHEMESTVD}.}\label{FIG020}
\end{center}
\end{figure}
\paragraph{Findings:}
Let us summarize the results from this test case of a simple Riemann problem.
The LF-solution is the jagged line in figure \ref{FIG018} and is neither accurate nor smooth.
LW-R, LW-MC and both WAF variations are accurate but suffer heavy oscillations at the discontinuities.
Since we will have to deal with shock waves in the supernova model, these schemes are not useful here.
% FORCE and GOD are of first order in accuracy only, which results in dispertion of the shock wave soon.
FORCE and GOD are of first order in accuracy only, which results in diffusion of the shock wave soon.
VL and ROE provide nice results for this test case, but still create oscillations.
Unfortunately it turned out that VL as given in chapter \ref{SCHEMESVL} breaks down in the spherical symmetric case and we will drop it here without any further discussion.
Concerning TVD schemes, the SUPERBEE limiter is slightly superior to the Van Leer limiter and for the sake of simplicity we will only use the former from now on.
For all schemes, the numbers $\epsilon^{(\infty)}_\rho$ are affected by round-offs at the discontinuities and are not very meaningful.
So there is not much difference between weighting $\mathbf f$ or $\mathbf w$.
We will only use the weighted flux from now on.
The conclusion is to pursue only three schemes:
\begin{itemize}
\item First-order centred scheme (FORCE)\\
representing the class of low accuracy, smooth solutions (CI)
\item Roe-scheme (ROE)\\
representing the class of high accuracy, oscillating solutions (CII)
\item Total variation diminishing weighted average flux scheme with SUPERBEE limiter (TVD-F-SA)\\
% representing the class of highly sophisticated solutions (CIII)
representing the class of high accuracy solutions with repressed oscillations (CIII)
\end{itemize}
Figure \ref{FIG020} demonstrates the high quality of TVD-F-SA, which will be favoured among the three schemes.
\subsection{Courant-Friedrichs-Lewy condition}\label{SECCFL}
The physical propagation speed must never exceed the numerical propagation speed on the grid, thus $v_{\mathrm{max}}\le\frac{\Delta x}{\Delta t}$.
For a fixed grid size, the time step $\Delta t$ is calculated by $\Delta t=\operatorname{CFL}\cdot\frac{\Delta x}{v_{\max}}$ with $\operatorname{CFL}\le 1$.
We have chosen $\operatorname{CFL}=0.99$ in the previous test case.
That the $\operatorname{CFL}$-parameter is of importance can be seen by the example of the LW-MC scheme, where $\operatorname{CFL}\le 0.98$ causes significant deviations as illustrated in figure \ref{FIG035}.
For ROE, SW and TVD-F-SA there is little influence of the $\operatorname{CFL}$-parameter, except the duration of the simulation, of course.
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig035.cfl.LW-MC.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig036.cfl.eps}
\end{tabular}
\caption{Left: LW-MC fails for $\operatorname{CFL}\le 0.98$ (plot for $\operatorname{CFL}=0.90,0.98,0.99$); Right: ROE, SW and TVD-F-SA for various $0.1\le\operatorname{CFL}\le 0.99$ (12 plots in total). The results of these schemes hardly depend on the $\operatorname{CFL}$-parameter.}\label{FIG035}
\end{center}
\end{figure}
\cleardoublepage
\newpage
.
\newpage
\section{Supernova models}\label{SEC005}
\renewcommand{\abstractname}{}
\begin{abstract}
{\it
After all the groundwork in the previous chapters, it's time to start with the main objectives.
Simulations will be carried out for all phases of the supernova. They will use different initial conditions, implement heat conduction, viscosity, etc.
The simulation results will be contrasted with the theoretical predictions made earlier.
This chapter will also cover numerical issues related to spherically symmetry such as appropriate discretization and the effect of higher-order methods.
% Nevertheless a framework will be setup at first, which will hold for all investigations.
Nevertheless a framework will be setup at first.
}
\end{abstract}
\subsection{Setup and nomenclature}\label{SETUP}
In the theory of differential equations, problems often can be scaled arbitrarily.
For the present model it is therefore beneficial to fix some parameters a priori.
In all subsequent tests an ejected mass by the supernova of three sun masses, an explosion energy of $10^{44}\,\mathrm{Joule}$ and a density of one particle per square centimetre is taken.
It is furthermore assumed that heat conduction becomes relevant after about $t_\mathrm{cool}=21300\,\mathrm{years}$.
$t_\mathrm{cool}$ will also be the time limit of all simulations dealing with the Sedov-Taylor phase only.
All tests will be performed assuming the equation of state of an ideal monoatomic gas with adiabatic index $\gamma=\frac{5}{3}$ and a mean molecular weight $\mu=\frac{13}{21}$ following from a hydrogen/helium mass ration of 3:1 and zero metallicity ($Z=0$).
For simplicity, dimensionless quantities as introduced in chapter \ref{SNR} are used.
\begin{align}
[L] &= 1\,\mathrm{pc} = 3.086\cdot 10^{16}\mathrm{m} &\qquad \text{unit of length} \\
[T] &= 10000\,\mathrm{a} = 3.156\cdot 10^{11}\mathrm{s} &\qquad \text{unit of time} \\
[M] &= 1\,\mathrm{M_\odot} = 1.989\cdot 10^{30}\mathrm{kg} &\qquad \text{unit of mass}
% [E] &= 10^{44}\,\mathrm{J} &\qquad \text{energy}
\end{align}
Therefore we have a fixed set of parameters as follows:
\begin{center}
\begin{tabular}{r|c|l}
fixed parameter & value & description \\
\hline
$E_\mathrm{SN}$ & $10^{44}\,\mathrm{J}$ & energy released by SN \\
$M^\ast$ & $3$ & total mass of ejecta \\
$\rho_\mathrm{ISM}^\ast$ & $0.0153$ & mass density of ISM \\
$v_\mathrm{ISM}^\ast$ & $0$ & ISM at rest \\
$\gamma$ & $\frac{5}{3}$ & adiabatic index for monoatomic gas \\
$\mu$ & $\frac{13}{21}$ & mean molecular weight \\
$T_\mathrm{ISM}$ & $10000\,\mathrm{K}$ & temperature of ISM \\
$t_\mathrm{cool}^\ast$ & $2.13$ & end of Sedov-Taylor phase \\
& & time when heat conduction becomes relevant \\
\end{tabular}
\end{center}
Every simulation starts at time $t_0$.
This is when the shock front is at position $r_0$.
The initial conditions behind the shock wave at $r_0$ are denoted by $\rho_0(r)$, $v_0(r)$, etc., while the quantities in front of the shock wave are labelled by $\rho_{\mathrm{ISM}}(r)$, $v_{\mathrm{ISM}}(r)$, etc.
We distinguish between different types of initial conditions.
\begin{itemize}
\item Sedov initial conditions:\\
The Sedov solution of the blast as given in chapter \ref{SEDOVSOLUTION} serves as initial condition for the simulation.
We differentiate between the solution in a medium with constant ($\omega=0$) density (SED0) and the solution in a medium with exponentially falling density ($\omega=-1.5$) in radial direction (SED15).
In any case, the original density distribution of the medium before the blast is continuous.
This is a very restrictive condition and, with the given parameters from chapter \ref{SETUP}, the initial position of the shock wave $r_0^\ast$ and the initial time $t_0^\ast$ are already fully determined.
The Sedov initial conditions are useful for long-term simulations.
\item Riemann-like initial conditions:\\
For the Sedov initial conditions, the initial position of the shock wave $r_0^\ast$ is too far from the origin to observe early effects such as the reverse shock.
It is therefore necessary to give up the condition of an originally continuous density.
The best way to do so is to fix $r_0^\ast$, rescale the Sedov solution from above to the range between the origin and $r_0^\ast$ and allow a step in the density at $r_0^\ast$.
The scaling factor $\lambda$ is given by
\begin{eqnarray}
\lambda=\frac{r_0^{\ast(\mathrm{RIE})}}{r_0^{\ast(\mathrm{SED})}}
\end{eqnarray}
and the appropriate initial time, velocity, density and pressure for the shock front at $r_0^{\ast\mathrm{(RIE)}}$ in the case $\omega=0$ (RIE0) are
\begin{eqnarray}
t_0^{\ast\mathrm{(RIE)}}=\lambda t_0^{\ast\mathrm{(SED)}} &\qquad&
\rho_0^{\ast\mathrm{(RIE)}}(r_0^{\ast\mathrm{(RIE)}})=\lambda^{-3}\rho_0^{\ast\mathrm{(SED)}}(r_0^{\ast\mathrm{(SED)}})\\
v_0^{\ast\mathrm{(RIE)}}=v_0^{\ast\mathrm{(SED)}} &\qquad&
p_0^{\ast\mathrm{(RIE)}}(r_0^{\ast\mathrm{(RIE)}})=\lambda^{-3}p_0^{\ast\mathrm{(SED)}}(r_0^{\ast\mathrm{(SED)}})
\end{eqnarray}
In the upcoming simulations $r_0^{\mathrm{(RIE)}}=0.5\,\mathrm{pc}$ will be used.
\item Specific density profiles:\\
It is not necessarily required to start with a precalculated blast solution.
Another way is to initially fix the density profile.
If this profile is meaningful, the simulation will adjust the other parameters accordingly.
One interesting profile is illustrated in \cite{DOR002} on page 35, where the density is small at the origin ($r=0$) and at the shell ($r=r_0$), while the peak in between is larger by 6--7 magnitudes.
Here we will adapt the profile given in equation (\ref{EQU040}).
We additionally require
\begin{itemize}
\item that all the ejected mass $M^\ast=3$ is within $r_0$
\item that the initial velocity field $v_0(r)=\sqrt{\frac{2E}{M}}\frac{r}{r_0}$ for $r\le r_0$ is linear and $v(r)=0$ for $r>r_0$
\item that at the beginning, continuity is given at $r_0$, thus $\rho_0(r_0)=\rho_\mathrm{ISM}$
\end{itemize}
This gives the uniquely determined initial density profile (IDP)
\begin{eqnarray}
\rho_0^\ast(r^\ast)=\begin{cases}
10.7993\cdot e^{-4.4267\cdot r^\ast}\qquad & r^\ast\le r_0^\ast\\
0.01531\qquad & r^\ast>r_0^\ast
\end{cases}
\end{eqnarray}
with $r_0^\ast=1.4817$ and $t_0^\ast=0.0250$.
\end{itemize}
\subsection{Source term}
All techniques applied here belong to the category of so-called finite volume methods.
These have the convenient property that all types of conservation equations maybe locally treated correctly.
The entire solution domain of the model is split by a grid into small cells $\Delta V$.
The homogeneous PDE
\begin{eqnarray}
% \frac{\partial}{\partial t}\mathbf{u}(t,x)+\vec\nabla\cdot f(\mathbf{u}(t,x))=s(\mathbf{u}(t,x))
\frac{\partial}{\partial t}\mathbf{u}(t,x)+\vec\nabla\cdot f(\mathbf{u}(t,x))=0
\end{eqnarray}
after applying Gauss's divergence theorem reads
\begin{eqnarray}
\int_{\Delta V}\frac{\partial}{\partial t}\mathbf{u}\,dV+\int_{\partial\Delta V} f(\mathbf{u})\cdot{\mathbf{n}}\,dS=0
\end{eqnarray}
This is exactly the idea behind these methods, namely that the temporal changes {\it within} the finite volumes $\Delta V$ are calculated due to flows through the surfaces $\partial\Delta V$ {\it enclosing} the volumes.
For the spherically symmetric case, these cells are naturally spherical shells and parameterized radially in only one dimension.
In the implementation presented here, the boundaries between the cells are exactly the integer multiples of the spatial discretization $\Delta x$, i.e., the fluxes $\mathbf f_{k+\frac{1}{2}}$ are located exactly at the points $x=k\Delta x$.
The point $x_k$, representing the interior of the cell $[(k-1)\Delta x,k\Delta x]$, is located somewhere between $x_\mathrm{L}=(k-1)\Delta x$ and $x_\mathrm{R}=k\Delta x$.
Scalar quantities such as density $\rho_k$ and pressure $p_k$ are defined at $x_k$, as shown in figure \ref{FIG037}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=15cm]{SN.fig037.eps}
% \includegraphics[width=\textwidth]{SN.fig037.eps}
\caption{Setup -- the fluxes $\mathbf{f}$ are located on integer multiples of the spatial resolution $\Delta x$.
The arrows indicate direction and amount of the flow.
The scalar components of $\mathbf{u}$ are the values of the unknowns in the cell $[(k-1)\Delta x,k\Delta x]$ and are not associated with a particular position.
To calculate the fluxes $\mathbf f$ no $x_k$ need to be specified.
To calculate the source term $g(\mathbf u)$ a consistent choice for the $x_k$ has to be made.}\label{FIG037}
\end{center}
\end{figure}
We have seen in chapter \ref{EULERGEOMETRIES} and chapter \ref{PDENUMERICALLY} how the geometry of the problem can be represented by source terms $s(\mathbf u)$.
In the case of spherical symmetry, the calculation splits into two steps, which are
\begin{enumerate}
\item Use a suitable differential scheme to estimate the fluxes $\mathbf{f}$ between the cells.
The general form for forwarding the scalars $\mathbf{u}$ within the cells always have the form
\begin{eqnarray}\label{EQU038}
\mathbf u_k^{(n+1)}=\mathbf u^{(n)}_k-\frac{\Delta t}{\Delta x}\left(\mathbf f_{k+\frac{1}{2}}-\mathbf f_{k-\frac{1}{2}}\right)
\end{eqnarray}
The exact position of $x_k$ is irrelevant for the calculation of the fluxes $\mathbf f_{k-\frac{1}{2}}$ and $\mathbf f_{k+\frac{1}{2}}$.
\item Forward the resultant ODE with suitable methods such as Runge--Kutta
\begin{eqnarray}\label{EQU039}
\frac{d\mathbf u^{(n+1)}}{dt}=g(\mathbf u^{(n+1)})=-\frac{2}{x}s(\mathbf u^{(n+1)})
\end{eqnarray}
Equation (\ref{EQU039}) requires to define the positions of the $x_k$.
In the planar case the obvious choice is
\begin{eqnarray}
x_k^{(M)}=(k-\frac{1}{2})\Delta x
\end{eqnarray}
but for the spherically symmetric case the center of mass
\begin{eqnarray}
x_k^{(C)}=\sqrt[3]{\frac{x_\mathrm{L}^3+x_\mathrm{R}^3}{2}}
\end{eqnarray}
with $x_\mathrm{L}=(k-1)\Delta x$ and $x_\mathrm{R}=k\Delta x$ seems to be more consistent.
\end{enumerate}
\subsubsection{Runge--Kutta}
To evaluate the effects of higher-order methods for the ODE caused by the source term, a first order simulation, a third order Runge--Kutta simulation and a fourth order Runge--Kutta simulation are compared with each other.
Setup parameters are given in table \ref{TAB004}.
% \begin{table}[H]
\begin{center}
\begin{tabular}{lclcl}
first order & & Runge--Kutta third order & & Runge--Kutta fourth order \\
\hline
$k_1 = g(\mathbf u)$ & &
$k_1 = g(\mathbf u)$ & &
$k_1 = g(\mathbf u)$ \\
& &
$k_2 = g(\mathbf u+\frac{\Delta t}{2}k_1)$ & &
$k_2 = g(\mathbf u+\frac{\Delta t}{2}k_1)$ \\
& &
$k_3 = g(\mathbf u-\Delta t\,k_1+2\Delta t\,k_2)$ & &
$k_3 = g(\mathbf u+\frac{\Delta t}{2}k_2))$ \\
& &
& &
$k_4 = g(\mathbf u+\frac{\Delta t}{2}k_3))$ \\
$\mathbf u \to \mathbf u+\Delta t\,k_1$ & &
$\mathbf u \to \mathbf u+ \frac{\Delta t}{6}(k_1+4k_2+k_3)$ & &
$\mathbf u \to \mathbf u+ \frac{\Delta t}{6}(k_1+2k_2+2k_3+k_4)$ \\
\end{tabular}
\begin{eqnarray}
g(\mathbf u_k) = -\frac{2}{x_k^{(C)}}s(\mathbf u_k)
\end{eqnarray}
\end{center}
% \caption{First, third and fourth order method to forward an ODE.}\label{TAB005}
% \end{table}
\begin{table}[H]
\begin{center}
\begin{tabular}{r|l}
Parameter & Value \\
\hline
initial condition & SED0 \\
scheme & ROE, TVD-F-SA \\
N & 3000 \\
$dx^\ast$ & 0.01 \\
$t_0^\ast$ & 0.0295 \\
$t_\mathrm{final}^\ast $ & 2.1295 \\
support points & center of mass (C)
\end{tabular}
\caption{Setup for evaluating the method to forward the source term.}\label{TAB004}
\end{center}
\end{table}
All tests performed in the framework of this thesis have shown that the influence of the method to forward the source term is negligible compared to the effects of the differential schemes to forward the homogeneous PDE.
% From now on we will use Runge--Kutta knowing that we are on the safe side with this choice.
\subsubsection{Support points}
We compare two identical simulations except for the support points $x_k$ which are located at
\begin{eqnarray}
x_k &=&
\begin{cases}
(k-\frac{1}{2})\Delta x\qquad & \text{case (M)} \\
\sqrt[3]{\frac{(k-1)^3+k^3}{2}}\Delta x=(k-\frac{1}{2}+\frac{1}{4k}+\frac{1}{8k^2}+\frac{5}{24k^3}+\dots)\Delta x\qquad & \text{case (C)}
\end{cases}
\end{eqnarray}
\begin{center}
\begin{tabular}{r|l}
Parameter & Value \\
\hline
initial condition & SED0 \\
scheme & TVD-F-SA \\
N & 3000 \\
$dx^\ast$ & 0.01 \\
$t_0^\ast$ & 0.0295 \\
$t_\mathrm{final}^\ast $ & 2.1295 \\
ODE method & Runge--Kutta fourth order
\end{tabular}
\end{center}
Again, all tests performed have shown that the effects of this parameterisation are negligible.
In this particular case, the maximum deviation at $t^\ast_{\mathrm{final}}$ is
\begin{eqnarray}
\frac{|\rho^{\ast(M)}-\rho^{\ast(C)}|}{4\rho_{\mathrm{ISM}}^\ast}<0.0001
\end{eqnarray}
We are going to use the method where the cell's center of mass is taken for $x_k$, since this approach is probably more consistent.
% Should: $r_{\mathrm{final}}^\ast=19.95$, IST: $r_{\mathrm{final}}^\ast=20.02$ (beide Fälle), $\frac{|\Delta\rho^\ast|}{4\rho_{\mathrm{ISM}}^\ast}<0.0001$\\
\newpage
\subsection{Initial conditions}
We will now compare the initial conditions SED0, SED15, RIE0 and IDP for the schemes FORCE, ROE and TVD-F-SA.
The parameters used are those introduced in chapter \ref{SETUP}.
The ODE is forwarded by third order Runge--Kutta with support points $x_k^{(C)}$.
The results at $t^\ast_{\mathrm{final}}=t^\ast_{\mathrm{cool}}=2.13$ are compared.
The durations of the individual simulations are consequently $t^\ast_{\mathrm{final}}-t^\ast_{0}$.
\begin{center}
\begin{tabular}{r|l|l|l|l}
& SED0 & SED15 & RIE0 & IDP \\
\hline
N & 3000 & 3000 & 3000 & 3000 \\
$dx^\ast$ & 0.01 & 0.01 & 0.01 & 0.01 \\
$t_0^\ast$ & 0.0295 & 0.0245 & 0.0041 & 0.0250 \\
$t_\mathrm{final}^\ast $ & 2.13 & 2.13 & 2.13 & 2.13 \\
$r_0^\ast$ & 3.605 & 2.865 & 0.495 & 1.485 \\
$\rho_0^\ast(r_0^\ast)$ & 0.06152 & 0.06171 & 20.3671 & 0.01531 \\
$v_0^\ast(r_0^\ast)$ & 36.6188 & 50.0988 & 35.8636 & 58.9408 \\
$p_0^\ast(r_0^\ast)$ & 27.4101 & 51.4862 & 9461.23 & 40.2415 \\
\end{tabular}
\end{center}
It is now time to use the theoretical predictions from chapter \ref{SNR}.
The density of the shock wave is expected to be $\frac{\gamma+1}{\gamma-1}\rho_\mathrm{ISM}=4\,\rho_\mathrm{ISM}$.
According to equation (\ref{EQU022}) its position at $t^\ast=t_\mathrm{cool}^\ast=2.13$ should be at about $r^\ast=19.95$.
The simulations discussed above produce the following results:
\begin{center}
\begin{tabular}{r|l|l|l|l}
$r^\ast$ & SED0 & SED15 & RIE0 & IDP \\
\hline
FORCE & 20.09 & 19.98 & 20.54 & 16.40 \\
ROE & 19.99 & 20.01 & 20.28 & 16.40 \\
TVD-F-SA & 20.04 & 20.06 & 20.12 & 16.57
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{r|l|l|l|l}
$\rho^\ast/4\,\rho_\mathrm{ISM}^\ast$ & SED0 & SED15 & RIE0 & IDP \\
\hline
FORCE & 0.699 & 0.885 & 0.828 & 0.923 \\
ROE & 0.948 & 0.956 & 0.946 & 0.944 \\
TVD-F-SA & 0.977 & 0.978 & 0.976 & 0.971 \\
\end{tabular}
\end{center}
Figures \ref{FIG050} and \ref{FIG060} illustrate the findings in a compact way.
For all four initial conditions (SED0, SED15, RIE0, IDP), a simulation with each scheme (FORCE, ROE, TVD-F-SA) was performed.
The curves are overlaid for comparison (9 curves for SED0, SED15, RIE0 in figure \ref{FIG050}, 3 curves for IDP in figure \ref{FIG060}).
For each quantity (density, velocity, pressure, internal energy), the position of the shock wave was enlarged to resolve the individual simulation results.
The previously discussed properties of the various differential schemes can be observed again.
In case of RIE0 and IDP initial conditions we find artifacts, which will be the topic of the next chapter.
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig050.rho.eps} &
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig052.v.eps} \\
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig054.p.eps} &
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig056.e.eps} \\
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig051.rho.eps} &
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig053.v.eps} \\
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig055.p.eps} &
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig057.e.eps}
\end{tabular}
\caption{Numerical solutions at $t^\ast=2.13$ for SED0, SED15 and RIE0 initial conditions for FORCE, ROE and TVD-F-SA schemes (9 curves).
The four bottom panels are magnifications of the four panels in the upper half of the figure.
The bump in the RIE0 initial condition case is due to a reverse shock and will be dealt with in the next chapter.}\label{FIG050}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig060.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig061.v.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig062.p.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig063.e.eps} \\
\end{tabular}
\caption{Numerical solutions at $t^\ast=2.13$ for the IPD initial condition.}\label{FIG060}
\end{center}
\end{figure}
The conclusions are:
\begin{itemize}
\item As an effect of diffusion, in all simulations the shock wave propagates slightly too fast.
The IDP initial conditions are apparently not applicable for the theoretical model.
\item The TVD-F-SA scheme gives best results for shape and amplitude of the shock wave.
The FORCE scheme is inaccurate and the shock wave dissolves considerably.
\end{itemize}
It should be mentioned that the low densities close to the origin have to be treated carefully in any scenario, even though they have hardly any impact on the system beyond $r_0$.
All simulations described above enforce a minimum density of $\rho_\mathrm{min}=10^{-6}\,\rho_\mathrm{ISM}$.
It has been verified that this modification has no effect on the propagation of the shock wave.
Also, the reverse shock has to be suppressed before it reaches the origin to avoid ambiguities.
This can be done, for example, by not allowing negative velocities close to the origin.
\subsection{Reverse shock}
To observe a significant reverse shock, a scheme of high order is required.
We will therefore use only the TVD-F-SA scheme with third order Runge--Kutta for the inhomogeneous term, and the centers of mass as support points of the finite volumes.
For SED0 and SED15 the initial position of the shock wave $r_0^\ast$ is too far from the origin to observe the reverse shock.
For RIE0 the initial position of the shock wave is at $r_0=0.495\,\mathrm{pc}$ and the reverse shock can be observed.
Figure \ref{FIG070} shows density and velocity between $t^\ast=0.2$ and $t^\ast=2.0$.
The jagged graphs of the velocities reveal that the system becomes chaotic in the domain of the reverse shock after a while.
Fortunately the propagating shock wave is in no causal dependence on this.
\begin{figure}[ht]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig070.rho.i3.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig072.v.i3.eps} \\
\end{tabular}
\caption{Reverse shock at a number of equally spaced points in time between $t^\ast=0.2$ and $t^\ast=2.0$ for RIE0 initial conditions.}\label{FIG070}
\end{center}
\end{figure}
The existence of the reverse shock is an established fact.
However, with the resolution used here, we cannot mirror all facets of a supernova exactly.
In particular, the simulation of the reverse shock should be seen as a proof of existence rather than an exact copy of reality.
We have seen that the IDP initial conditions describe the supernova model only with reservations.
For studying the reverse shock, however, this setup exhibits nice results, as shown in figure \ref{FIG074}.
A variation of the setup by using equation (\ref{EQU040}) is illustrated in figure \ref{FIG076}.
Here the reverse shock accelerates until the simulation breaks down numerically in a region around the origin at about $t^\ast=0.48$.
\begin{figure}[ht]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig074.rho.i5.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig075.v.i5.eps} \\
\end{tabular}
\caption{Reverse shock at a number of equally spaced points in time for $t^\ast=0.10, 0.15,\dots, 0.50$ for IDP initial conditions.}\label{FIG074}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig076.rho.i6.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig077.v.i6.eps} \\
\end{tabular}
\caption{Reverse shock at a number of equally spaced points in time for $t^\ast=0.1, 0.2,\dots, 1.0$ for a variation of IDP initial conditions using equation (\ref{EQU040}). The reverse shock hits the origin at about $t^\ast=0.48$.}\label{FIG076}
\end{center}
\end{figure}
\subsection{Heat conduction}
% Admittedly the analysis of the reverse shock is a very vague undertaking.
% Heat conduction in comparison is easier to handle.
% The effects of cooling are only relevant in regions far from the origin and we are not facing the collision of waves.
Admittedly the analysis of the reverse shock has been a rather qualitative discussion.
By contrast, heat conduction should be easier to handle.
The effects of cooling are only relevant in regions far from the origin and it is not necessary to account for the collision of waves.
Cooling is affected by a heat conduction term in the energy equation as already introduced in equation (\ref{EQU008}).
As in the case of the source term, the homogeneous Euler equations are treated first.
After that, using the new values for $\mathbf u$, the ODE
\begin{eqnarray}
\frac{d\epsilon^{(n+1)}}{dt}&=&-\left(\frac{\rho^{(n+1)}}{\mu m_\text{H}}\right)^2\Lambda(T(\rho^{(n+1)},p^{(n+1)})) \\
T(\rho^{(n+1)},p^{(n+1)}) &=& \frac{\mu m_\text{H}}{k_B}\cdot\frac{p^{(n+1)}}{\rho^{(n+1)}}
\end{eqnarray}
has to be forwarded in time.
% As can be seen in figure \ref{FIG001}, the cooling function $\Lambda(T)$ covers for relevant $T$ a range of two orders of magnitude.
As can be seen in figure \ref{FIG001}, in the relevant temperature range the cooling function $\Lambda(T)$ varies by two orders of magnitude.
The consequence is that at approximately $300000\,\mathrm{K}$ large gradients in the internal energy $\epsilon$ emerge which result in non-physical densities behind the shock wave.
This problem can be attacked by a finer grid.
However, this is not feasible since calculation time increases at least quadratically with grid size.
The alternative way to proceed (at least within this framework) is to limit the gradient.
\begin{eqnarray}
\frac{d\epsilon}{dt}=-\min\left(\beta\epsilon,n^2\Lambda\right)
\end{eqnarray}
This is justified by the fact that the cooling maximum at $\log_{10}\mathrm{T}/\mathrm{K}=5.5$ is, among other things, caused by the presence of the elements $\mathrm{C}$ and $\mathrm{O}$.
Since we have assumed metallicity $Z=0$, these two elements do not appear in the present model and the appropriate cooling should be lower in this range.
Using the SED0 initial conditions from chapter \ref{SETUP} on a grid with resolution $\Delta x^\ast=0.01$ for all values of $\beta$ in the interval $[0.25,0.75]$ the numerical problem vanishes;
otherwise $\beta<0.25$ causes the shock front to diffuse and $\beta>0.75$ cannot resolve the wave correctly.
We will proceed with $\beta=0.5$ knowing that results for any other value $0.25\le\beta\le 0.75$ do not differ substantially.
The simulation output is not very surprising.
The propagation speed of the shock front slows down at a faster rate than in the solution without heat conduction (see table \ref{TAB003}) and the density within the wave piles up.
Figure \ref{FIG080} shows the solution for TVD-F-SA for $t^\ast=5,10,\dots,50$.
The width of the shock front is diminished significantly and as a consequence the density $\rho_\mathrm{cool}^\ast$ is up to 20 times larger than the density maximum without cooling ($\rho^\ast\approx 0.06$).
It is interesting to observe the temperature minimum behind the front which is caused by the strong cooling at around $300000\,\mathrm{K}$.
Finally we compare the position of the shock front with the predictions made in chapter \ref{SNR}, where we deduced the exponential correlation between the actual simulation time $t^\ast$ and the position $r^\ast$ of the shock front
\begin{eqnarray}
{r^\ast} &\propto& {t^\ast}^\alpha
\end{eqnarray}
with $\alpha=\frac{2}{5}$ in the case without cooling and $\alpha=\frac{2}{7}$ in the case when cooling is applied.
These results are summarized in table \ref{TAB003}.
In both cases, the shock front in the simulation moves slightly faster than expected.
This can be explained by a dissipative mechanism in Godunov-type schemes \cite{XU001}.
\begin{table}[H]
\begin{center}
\begin{tabular}{c|ccc|ccc}
$t^\ast$ & simulation & prediction & error & simulation & prediction & error \\
cooling & & NO & & & YES & \\
\hline
5 & 28.28 & 28.06 & +0.77\% & 26.08 & 25.52 & +2.19\% \\
10 & 37.31 & 37.03 & +0.75\% & 32.42 & 31.11 & +4.21\% \\
15 & 43.92 & 43.55 & +0.85\% & 36.89 & 34.93 & +5.60\% \\
20 & 49.32 & 48.86 & +0.94\% & 40.36 & 37.93 & +6.42\% \\
25 & 53.99 & 53.42 & +1.06\% & 43.12 & 40.42 & +6.67\% \\
30 & 58.14 & 57.47 & +1.17\% & 45.42 & 42.58 & +6.66\% \\
35 & 61.92 & 61.12 & +1.31\% & 47.42 & 44.50 & +6.56\% \\
40 & 65.40 & 64.47 & +1.44\% & 49.22 & 46.23 & +6.46\% \\
45 & 68.65 & 67.58 & +1.58\% & 50.80 & 47.81 & +6.25\% \\
50 & 71.70 & 70.49 & +1.71\% & 52.23 & 49.28 & +6.00\%
\end{tabular}
\caption{Position of the shock front at $t^\ast$ with and without cooling -- simulations versus predictions based on astrophysical considerations done in chapter \ref{SNR}.
The TVD-F-SA scheme with SED0 initial conditions is used on a grid with resolution $\Delta x^\ast=0.01$.
}\label{TAB003}
\end{center}
\end{table}
\begin{figure}[H]
\begin{center}
\includegraphics[width=13.8cm,height=8.5cm]{SN.fig080.eps}
\caption{Density $\rho^\ast$ of the shock wave within the first $0.5\,\mathrm{Myrs}$ when cooling is applied.
TVD-F-SA scheme with SED0 initial conditions on a grid with resolution $\Delta x^\ast=0.01$.
}\label{FIG080}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=13.8cm,height=8.5cm]{SN.fig081.eps}
\caption{Temperature of the shock wave within the first $0.5\,\mathrm{Myrs}$ when cooling is applied.
TVD-F-SA scheme with SED0 initial conditions on a grid with resolution $\Delta x^\ast=0.01$.
}\label{FIG080a}
\end{center}
\end{figure}
\subsection{Artificial viscosity}\label{SECARTIFICALVISCOSITY}
Simulations of shock waves occasionally require the inclusion of artificial viscosity for damping the oscillations introduced by the numerical scheme.
One sophisticated way is given in \cite{CIO001} by adapting the second Euler equation (momentum equation) using equation (\ref{EQU042}) as basis.
Unfortunately this approach cannot be implemented straightforwardly for most schemes including ROE and TVD, but needs extensive modifications.
We follow a simpler idea suggested in 1967 by Richtmyer \& Morton \cite{RIE001} (p313).
Here a viscosity-like term $q$ is added to the pressure.
This gives the modified Euler equations
\begin{eqnarray}
\partial_t\rho &=& -\vec\nabla(\rho\mathbf v) \\
\partial_t(\rho\mathbf v) &=& -\vec\nabla((\rho\mathbf v)\mathbf v+(p+q)\mathbf I) \\
\partial_t\epsilon &=& -\vec\nabla(\epsilon\mathbf v+(p+q)\mathbf v) \\
q &=&
\begin{cases}
l^2\rho\left(\frac{\partial v}{\partial r}\right)^2 & \frac{\partial v}{\partial r}<0,\qquad [l]=L \\
0 & \frac{\partial v}{\partial r}\ge 0 \\
\end{cases}
\end{eqnarray}
where $l$ is a parameter of dimension length, which has to be chosen adequately.
This approach still needs considerable adaption in the weighting of the quantities for the complex TVD-schemes, but works fine for most other schemes discussed in this thesis.
For most of the initial conditions considered here, the initial value of $\frac{\partial v}{\partial r}$ at the shock front is infinite.
With decreasing $\Delta x$ the numerical estimate becomes very large.
Therefore, just like for the cooling term $n^2\Lambda$, $q$ also needs to be limited
\begin{eqnarray}
q_\mathrm{lim} &=&\min(q,\beta p)
\end{eqnarray}
Trial simulations suggest $\beta=0.5$ as reasonable value when using $l^2=0.5$.
Figure \ref{FIG082} shows the effect of adding the viscosity term $q$ for the ROE-scheme.
% The improvement of numerical stability has to be paid by a diffusion of the wave front.
The improvement of numerical stability is obtained at the price of blurring the wave front.
Artificial viscosity is indisputably a powerful tool in fluid dynamics.
Yet for the supernova simulation it seems to be advantageous to invest into the complexity of the scheme rather than into this extension.
\begin{figure}[H]
\begin{center}
\includegraphics[width=13.8cm,height=8.5cm]{SN.fig082.eps}
\caption{Artificial viscosity: numerical stability has to be paid for by a diffusive/dissolving shock wave.
TVD-F-SA scheme with SED0 initial conditions on a grid with resolution $\Delta x^\ast=0.01$.
}\label{FIG082}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=13.8cm,height=8.5cm]{SN.fig083.eps}
\caption{Artificial viscosity: $q$ needs to be limited where $\frac{\partial v}{\partial r}$ becomes large.
TVD-F-SA scheme with SED0 initial conditions on a grid with resolution $\Delta x^\ast=0.01$.
}\label{FIG083}
\end{center}
\end{figure}
\newpage
\section{A roadmap for the future}\label{SEC006}
\renewcommand{\abstractname}{}
\begin{abstract}
{\it
This last chapter starts with a brief summary of what has been achieved.
What follows is a variety of themes that promise a wide range of possible activities for the future.
}
\end{abstract}
\paragraph{Summary:}
Let us summarize what has happened so far.
First of all we started by comparing various numerical methods for solving nonlinear systems of differential equations.
We realized that higher-order methods improve the accuracy of the solution, but at the same time cause problems such as oscillations.
To counteract this behaviour we discussed a group of differential schemes based on the principle of {\it total variation diminishing}.
The so-called {\it flux limiter functions} with {\it flow parameter} arguments are used to smooth out numerical artifacts in the solution.
To test these differential schemes we set up an astrophysical model case, namely the propagation of a supernova.
In its simplest form it involves the spherically symmetric solution of the homogeneous Euler equations after a point blast.
We confirmed our simulation results with the idealized model solutions and simultaneously enhanced the methodology.
Special attention was given to the source term, an inhomogeneous term in the Euler equations which occurs due to the spherical symmetry.
The next step consequently was the extension of the model.
We investigated different initial conditions, examined the reverse shock after the free expansion and simulated heat conduction by inserting a cooling term.
However, we realized that with the given relatively simple model the limits of what is possible are reached very quickly.
\paragraph{Perspectives:}
% I have to admit that all of above was just a humble beginning for a complete description of a supernova.
I see this thesis as only a humble beginning of a complete description of a supernova.
How do we go on from here?
A first step on the roadmap is certainly the relaxation of spherical symmetry.
This is necessary to simulate the instabilities occurring near the shock front.
It will furthermore allow us to deal with inhomogeneities in the ISM, such as changing density due to e.g.\ molecular clouds.
Of particular interest should be the collision of two supernova fronts.
The resulting mass concentrations are known to be a prerequisite for new star formation.
Another important issue which has been neglected so far is the initial phase just after the blast.
What all simulations within this thesis have in common is that they start in the late stage of the free expansion phase.
A comprehensive model would also need to consider the very early stage, where particle interactions require an alternative approach to the Euler equations.
With such modifications it should be possible to introduce metallicity into the model as well.
Also on the methodological side there are many interesting questions.
% It might be advantageous to use a Lagrangian description instead of the Euler approach.
It might be advantageous to use a Lagrangian description by tracking the motion of the fluid particles, rather than the present Euler description where fluid properties are functions of space and time.
Particularly a smoothed particle hydrodynamics (SPH) model seems to be promising.
Large simulations with more features use much finer grids, which are only feasible with extensive performance optimization.
In fact, this is an area of research by itself.
\paragraph{Final words:}
I would like to conclude this work with a comparison.
For me, the process of studying the subject of supernovae has been similar to that of learning a piano piece.
In both cases, one spends considerable amount of time until there is an initial understanding and satisfaction.
% What has been accomplished is then allowed to rest, so it can be taken up again after a while to look at it from a different angle.
The work is then put away for a while, so it can be taken up again after some time to look at it from a different angle.
In the process of writing this thesis this happened to me at least three times, and each time I gained a new perspective.
This was extremely beneficial in improving my understanding of the subject and consequently the quality of the result.
% I hope to be able to pick up some of the open issues while working on a prolongation of this thesis.
% my astronomy thesis, which I see as a continuation of this work.
I hope to be able to pick up some of the open issues during my future academic pursuits.
\paragraph{Addendum:}
A sophisticated perspective for the future using the techniques outlined in this thesis is simulating the collision of two supernova shock fronts.
This causes the agglomeration of density as required for new star formation.
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig094.rho.eps} &
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig095.v.eps} \\
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig096.p.eps} &
\includegraphics[width=7.8cm,height=4.8cm]{SN.fig097.e.eps} \\
\end{tabular}
\caption{The collision of two supernova shock fronts at $t^\ast=-0.25$, $t^\ast=0.05$ and $t^\ast=0.35$ where $t^\ast=0$ is the time of the collision.
Page 60: two shock waves approaching each other at $t^\ast=-0.25$.
Page 61, upper half: shortly after the collision at $t^\ast=0.05$. The density piles up, pressure and internal energy decrease.
Page 61, lower half: At $t^\ast=0.35$ the pressure in the interior of the supernova remnants has formed new shock waves, which are pushed back by the density at the point of collision.}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7.8cm,height=5cm]{SN.fig098.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig099.v.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig100.p.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig101.e.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig090.rho.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig091.v.eps} \\
\includegraphics[width=7.8cm,height=5cm]{SN.fig092.p.eps} &
\includegraphics[width=7.8cm,height=5cm]{SN.fig093.e.eps} \\
\end{tabular}
% \caption{}
\end{center}
\end{figure}
\newpage
.
\newpage
\appendix
\section{Program documentation}\label{SEC007}
The source code for the supernova simulations is written in $\mathrm{FORTRAN}$.
Program parameters are assigned in C-style via the command line.
The following table gives an overview of the input parameters.
\subsection{Parameters}
% \begin{table}
\begin{center}
\begin{tabular}{r|l}
Parameter=default & description \\
\hline
scheme=13 & defines the scheme to be used \\
& = 1 \dots exact Riemann solver \\
& = 2 \dots Lax--Friedrichs (LF) \\
& = 3 \dots Lax--Wendroff Richtmeyr (LW-R) \\
& = 4 \dots Godunov (GOD) \\
& = 5 \dots First-Order Centred Scheme (FORCE) \\
& = 6 \dots van Leer \\
& = 7 \dots Steger--Warming (SW) \\
& = 8 \dots Roe (ROE) \\
& = 9 \dots MacCormack (LW-MC) \\
& = 11 \dots WAF, flux averaged \\
& = 12 \dots WAF, $\mathbf w$ averaged \\
& = 13 \dots WAF TVD, flux averaged \\
& = 14 \dots WAF TVD, $\mathbf w$ averaged \\
& = 15 \dots MUSCL \\
\hline
init=4 & defines initial condition \\
& = 1 \dots Sedov, constant density (SED0) \\
& = 2 \dots Sedov, exp.decreasing density (SED15) \\
& = 3 \dots Riemann-like, constant density (RIE0) \\
& = 5 \dots decreasing cont. density profile (IDP) \\
& = 6 \dots equation (\ref{EQU040}) \\
& = 11 \dots the test case from chapter \ref{SEC004} \\
\hline
n=3000 & number of spatial points \\
\hline
dx=0.01 & length interval $\Delta x$ \\
\hline
t=2.1 & total run time in $10^4$ years \\
\hline
sym=2 & defines geometry \\
& = 0 \dots planar \\
& = 1 \dots cylindrical symmetry \\
& = 2 \dots spherical symmetry \\
\hline
limiter=1 & limiter function \\
& = 0 \dots no limiter function \\
& = 1 \dots SUPERBEE limiter \\
& = 2 \dots Van Leer limiter \\
& = 3 \dots Van Albada limiter \\
\hline
m=3.0 & mass of supernova in $M_\odot$ \\
\hline
\end{tabular}
\end{center}
% \caption{List of parameters, part one}
% \end{table}
\begin{center}
\begin{tabular}{r|l}
Parameter=default & description \\
\hline
CFL=0.99 & Courant-Friedrichs-Lewy number \\
\hline
RK=4 & Runge--Kutta for source term \\
& = 0, 1 \dots no Runge--Kutta, first order \\
& = 3 \dots Runge--Kutta third order \\
& = 4 \dots Runge--Kutta fourth order \\
\hline
cool=1 & cooling \\
& = 0 \dots cooling off \\
& = 1 \dots cooling on \\
\hline
v=3 & verbose parameter, defines amount of output \\
& = 1 \dots final result only \\
& = 3 \dots first and last state \\
& = 4 \dots states in predefined time-intervalls \\
& = 5 \dots every state \\
\hline
out=1.0 & time interval between output of intermediate results \\
& = 0.0 \dots no output of intermediate results \\
\hline
\end{tabular}
\end{center}
% \hline
% smear=0 & smearing discontinuity \\
% & gives width of smearing \\
% \hline
% CMB=0 & Cioffi, McKee \& Bertschinger \\
% & = 0 \dots OFF \\
% & = 1 \dots ON \\
% \hline
% \newpage
\newpage
\section{Curriculum Vitae}
\subsection*{Personal data}\nonumber
\begin{tabular}{p{48mm}p{100mm}}
\hline
Name & Dipl.-Ing. Dr. Michael Kenn, BSc \\
& *August $1^\mathrm{st}$, 1969, Vienna \\
Spouse & Ritu, born in New Delhi \\
Children & Shonali (*1999) and Sanjay (*2001) \\
Languages & German, English \\
\end{tabular}
\subsection*{Education}\nonumber
\begin{tabular}{p{48mm}p{100mm}}
\hline
High School (1980--1987) & Silver and Bronze medallist in national and international Math Olympics \\
Mathematics (1987--1992) & Dipl.-Ing. in Technical Mathematics at the Vienna University of Technology,
emphasis on information and data processing, received scholarship for excellent academic results,
graduation exam in biomathematics and number theory, passed with honors \\
Doctorate (1992--1994) & Dissertation at the Institute of Inorganic Chemistry on mathematical methods in IR-spectroscopy,
doctoral exam in analysis and functional analysis \\
Project Management (1998) & Participation in the Philips high potential program,
including leadership \& coaching, personal development and character building training \\
Astronomy (2009--2012) & BSc in Astronomy at the Vienna University,
passed with honors, received merit scholarships \\
Space Studies (2011) & Student at the International Space University, emphasis on Space Physics \\
Physics (2013--present) & MSc in Physics at the Vienna University, emphasis on Computational Physics, Relativity \& Cosmology and Particle Physics,
received merit scholarship \\
\end{tabular}
\subsection*{Professional experience}\nonumber
\begin{tabular}{p{48mm}p{100mm}}
\hline
Philips (1994--2008) & Mathematician at Philips Speech Recognition Systems \\
& 1994--1997:
Speech recognition technologist with focus on language modeling \\
& 1998--2001:
Group leader, responsible for the development of new languages \\
& 2002--2008:
Project manager for analysis, development and improvement of speech recognition components in the department for technology and research
\end{tabular}
\subsection*{Jobs and activities}\nonumber
\begin{tabular}{p{40mm}p{108mm}}
\hline
1987, 1989 & Summer internship at IBM Austria \\
1990 & Tobacco priming in Ontario, Canada \\
1990--1991 & Supervisor of math tutorials for electric engineers \\
1991--1992 & Professional ski instructor training \\
1991--1998 & Station manager in Val d'Isere and St.\ Anton at Freizeit Aktiv \\
1992--1993 & Teaching assistant at the Institute of Inorganic Chemistry \\
1997--1998 & Military service, worked as an expert on spreading patterns of infections at the ABC-defence school of the Austrian military \\
2002 & Paternity leave, stay in Goa/India for triathlon training \\
2013--2014 & Technical assistant at the Center of Medical Statistics, Informatics, and Intelligent Systems (CeMSIIS) \\
\end{tabular}
\subsection*{Personal interests}\nonumber
\begin{tabular}{p{40mm}p{108mm}}
\hline
Sports & Start in 300+ national and international competitions in various sports and countries \\
Travelling & Bicycle trips in 38 countries so far \\
Science & Participation in international number theory projects \\
Music & Piano with focus on classical music
\end{tabular}
\subsection*{Contact}\nonumber
\begin{tabular}{p{40mm}p{108mm}}
\hline
Telephone & +43-1-714\,18\,84 \\
e-mail & \href{mailto:michael@kenn.at}{michael@kenn.at}
\end{tabular}
\newpage
\section{Zusammenfassung}
Die hier vorliegende Arbeit behandelt das numerische Lösen von nichtlinearen partiellen Diffe\-rentialgleichungssystemen.
Am Beispiel einer Supernovaexplosion werden verschiedene Lösungs\-ansätze mit astrophysikalischen Modellen verglichen und optimiert.
Die besten Verfahren werden danach verwendet, um einfache Fragestellungen wie beispielsweise die Kühlung oder die Ausbrei\-tungsgeschwindigkeit der Schockwelle für verschiedene Anfangsbedingungen zu studieren.
Zunächst wird ein hydrodynamisches Modell für die Ausbreitung der Reste nach einer Supernova\-explosion beschrieben.
Die zugrundeliegenden Gleichungen sind dabei die Euler--Gleichungen.
Diese werden um inhomogene Terme erweitert, die Kühlung, Viskosität und dergleichen beschreiben.
Durch die riesige Menge an freigesetzter Energie entsteht eine Schockwelle, die bei ihrer Ausbrei\-tung ins interstellare Medium mehrere Phasen durchläuft.
Diese Phasen werden astrophysikalisch beschrieben und mathematisch formuliert.
Im mathematisch-numerischen Teil wird zunächst kurz auf die Theorie der partiellen Differentialglei\-chungen eingegangen.
Dabei werden Konsistenz, Stabilität, Genauigkeit und Konvergenz erörtert.
Es werden exakte Lösungen sowohl für das Riemann-Problem also auch für das Sedov-Problem einer punktförmigen Explosion beschrieben.
Im Kernteil werden dann 13 verschiedene Differenzen\-schemen miteinander verglichen.
Es stellt sich heraus, dass mit steigender Genauigkeit der Trend zu Oszillationen in der Lösung zunimmt.
Als eine mögliche Methode, dieses Problem zu beheben, wird das Prinzip der totalen Variationsminimierung ({\it total variation diminishing}) eingeführt.
Dabei werden sogenannte {\it flux limiter}-Funktionen verwendet, deren Argumente {\it flow parameter} sind und die eine Glättung der Lösungen bewirken.
Es folgt eine Vielzahl an Supernovasimulationen.
Die homogenen Euler--Gleichungen werden mit verschiedenen Typen von Differenzenschemen gelöst, und es wird auf die verschiedenen Eigenschaften der Lösungen eingegangen.
Für die inhomogenen Terme werden Runge--Kutta Verfahren verschiedener Ordnung miteinander verglichen.
Ebenso werden verschiedene Positionen der Stützstellen für diese Verfahren getestet.
Da die Simulationen erst unmittelbar nach der Phase der freien Expansion starten, also erst ein paar hundert Jahre nach der Explosion, wird speziell den verschiedenen möglichen Anfangsbedingungen besondere Aufmerksamkeit geschenkt.
Was folgt sind Untersuchungen zum {\it reverse shock} unmittelbar nach der freien Expansion, zur Kühlung der Schockfront, zur Ausbrei\-tungsgeschwindigkeit und schließlich noch die Einführung einer künstlichen Viskosität als zusätz\-licher Parameter.
Die Arbeit schließt mit einer Zusammenfassung und einem Ausblick auf notwendige Modellerweiterungen, die in einer Folgearbeit behandelt werden könnten.
\newpage
.
\newpage
% \tableofcontents
\begin{thebibliography}{99}
\bibitem{BAS001}
P.~Basieux,
{\it Die Top Seven der mathematischen Vermutungen},
rororo, Reinbek bei Hamburg
September 2004
\bibitem{CIO001}
D.~F.~Cioffi, C.~F.~McKee, E.~Bertschinger,
{\it Dynamics of Radiative Supernova Remnants},
Astrophysical Journal, Part 1 (ISSN 0004-637X), Vol. 334, Nov. 1, 1988, p. 252--265
\bibitem{DOR004}
E.~A.~Dorfi, L.~O'C.~Drury,
{\it Simple Adaptive Grids for 1-D Initial Value Problems},
Journal of Computational Physics, Vol. 69, Issue 1, March 1987, p. 175--195
\bibitem{DOR003}
E.~A.~Dorfi, H.~Pikall, A.~Stökl, A.~Gautschy,
{\it Towards a more consistent discretization scheme for adaptive implicit RHD computations},
Computer Physics Communications 174 (2006) p. 771--782
\bibitem{DWA001}
V.~V.~Dwarkadas, R.~A.~Chevalier,
{\it Interaction of Type IA Supernovae with Their Surroundings}
The Astrophysical Journal, Vol. 497, Issue 2, p. 807--823., 1998,\\
\url{http://adsabs.harvard.edu/abs/1998ApJ...497..807D}
% , ApJ, 497, 807, 1998
\bibitem{FIE001}
G.~B.~Field, D.~W.~Goldsmith, H.~J.~Habing,
{\it Cosmic-Ray Heating of the Interstellar Gas},
Astrophysical Journal, Vol. 155, p. L149,
1969
\bibitem{MCK001}
C.~F. McKee, E.~C.~Ostriker,
{\it A theory of the interstellar medium - Three components regulated by supernova explosions in an inhomogeneous substrate},
Astrophysical Journal, Part 1, Vol. 218, Nov. 15, 1977, p. 148--169.\\
{\scriptsize
\url{http://adsbit.harvard.edu/cgi-bin/nph-iarticle_query?bibcode=1977ApJ...218..148M}
}
\bibitem{MCK002}
C.~F. McKee, E.~C.~Ostriker,
{\it Theory of Star Formation},
Annual Review of Astronomy and Astrophysics, Vol. 45, p. 565--687,
September 2007\\
\url{http://www.annualreviews.org/doi/abs/10.1146/annurev.astro.45.051806.110602}
\bibitem{LAN001}
L.~D.~Landau, E.~M.~Lifschitz,
{\it Fluid Mechanics},
Pergamon Press, Oxford--New York--Beijing--Frankfurt--Sao Paulo--Sydney--Tokyo--Toronto,
% 1959/1987
1959
\bibitem{LEV001}
R.~R.~LeVeque,
{\it Finite-Volume Methods for Hyperbolic Problems},
Cambridge Texts in Applied Mathematics,
2002
\bibitem{POZ001}
C.~Pozrikidis,
{\it Introduction to Theoretical and Computational Fluid Dynamics},
Oxford University Press, $2^{\mathrm{nd}}$ edition,
September 28, 2011
\bibitem{RIE001}
R.~D.~Richtmyer, K.~W.~Morton,
{\it Difference Methods for Initial-Value Problems},
Krieger Pub Co, $2^{\mathrm{nd}}$ edition, New York--London--Sydney,
June 1994
\bibitem{ROE001}
P.~L.~Roe,
{\it Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes},
Journal of Computational Physics, Vol. 43, Issue 2, October 1981, p. 357--372\\
\url{http://www.sciencedirect.com/science/article/pii/0021999181901285}
\bibitem{SCN001}
P.~Schneider,
{\it Extragalaktische Astronomie und Kosmologie},
Springer Spektrum, Bonn, korr. Nachdruck 2007
\bibitem{SED001}
L.~I.~Sedov,
{\it Examples of Gas Motion and Certain Hypotheses on the Mechanism of Stellar Outbursts},
Reviews of Modern Physics, Vol. 30, Issue 3, p. 1077--1079,
1958\\
\url{http://adsabs.harvard.edu/abs/1958RvMP...30.1077S}
\bibitem{SED002}
L.~I.~Sedov,
{\it Similarity and dimensional methods in mechanics},
Academic Press, New York,
$10^{\mathrm{th}}$ edition, 1993
\bibitem{SHA001}
M.~Shabouei, R.~Ebrahimi, Ku.~Mazaheri,
{\it Numerical Simulation of Converging Spherical Shock Problem},
Ninth International Congress of Fluid Dynamics \& Propulsion (ICFDP9), December 18--21, 2008, Alexandria, Egypt \\
\url{http://www.icfd11.org/ICFDP9/ICFDP9-EG-267.pdf}
\bibitem{SUT001}
R.~S.~Sutherland, M.~A.~Dopita,
{\it Cooling functions for low-density astrophysical plasmas},
Astrophysical Journal Supplement Series (ISSN 0067-0049), Vol. 88, No. 1, p. 253--327,
1993\\
\url{http://adsabs.harvard.edu/cgi-bin/bib_query?1993ApJS...88..253S}
\bibitem{TAY001}
G.~Taylor,
{\it The Formation of a Blast Wave by a Very Intense Explosion - I. Theoretical Discussion},
Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, Vol. 201, No. 1065. (Mar. 22, 1950), p. 159--174 \\
{\scriptsize
\url{http://www-rohan.sdsu.edu/~jmahaffy/courses/f09/math636/lectures/allometric_dim/taylor.blast.wave.I.pdf}
}
\bibitem{TAY002}
G.~Taylor,
{\it The Formation of a Blast Wave by a Very Intense Explosion. II. The Atomic Explosion of 1945},
Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, Vol. 201, No. 1065. (Mar. 22, 1950), p. 175--186 \\
{\scriptsize
\url{http://www-rohan.sdsu.edu/~jmahaffy/courses/f09/math636/lectures/allometric_dim/taylor.blast.wave.II.pdf}
}
\bibitem{TOR001}
E.~F.~Toro,
{\it Riemann Solvers and Numerical Methods for Fluid Dynamic},
Springer, Berlin--Heidelberg, $3^{\mathrm{rd}}$ edition (April 27, 2009)
\bibitem{TOR002}
E.~F.~Toro,
{\it The Weighted Average Flux Method Applied to the Euler Equations},
Phil. Trans. R. Soc. Lond. A 1992 341, doi:10.1098/rsta.1992.0113, published 15 December 1992
\url{http://rsta.royalsocietypublishing.org/content/341/1662/499.full.pdf}
\bibitem{TRA001}
J.~A.~ Tragenstein,
{\it Numerical Solution of Hyperbolic Partial Differential Equations},
Cambridge University Press, 2009
\bibitem{TRA002}
J.~A.~ Tragenstein,
{\it Numerical Solution of Elliptic and Parabolic Partial Differential Equations},
Cambridge University Press, 2013
\bibitem{XU001}
Kun Xu, Zuowu Li, {\it Dissipative mechanism in Godunov-type schemes},
International Journal for Numerical Methods in Fluids,
Vol. 37, p. 1--22 (doi:10.1002/fld.160), 2001
\section*{Non-refereed References}\nonumber
\bibitem{DOR001}
E.~A.~Dorfi, {\it Vorlesungsskript Astrophysik I}, Universität Wien, 2011
\bibitem{DOR002}
E.~A.~Dorfi, {\it Vorlesungsskript Astrophysik II}, Universität Wien, 2012
\bibitem{DWA002}
V.~V.~Dwarkadas, {\it The Interaction of Type Ia Supernovae with their Surroundings: The Exponential profile in 2D}, 2006\\
\url{http://cds.cern.ch/record/450217/files/0008076.pdf}
\bibitem{KAM001}
J.~R.~Kamm, {\it Evaluation of the Sedov-von Neumann-Taylor Blast Wave Solution}, 2000\\
\url{http://cococubed.asu.edu/papers/kamm_2000.pdf}
\bibitem{KAM002}
J.~R.~Kamm, F.~X.~Timmes, {\it On Efficient Generation of Numerically Robust Sedov Solutions} \\
\url{http://cococubed.asu.edu/papers/la-ur-07-2849.pdf}
\bibitem{NEU003}
M.~Neumann, {\it Computational Physics Vorlesungsskript - Partielle Differentialgleichungen}\\
\url{http://www.exp.univie.ac.at/cp1/cp3.pdf}
\bibitem{NIS001}
H.~Nishikawa, {\it A Comparison of Numerical Flux Formulas for the Euler Equations}, 1998\\
\url{http://www.cfdnotes.com/cfdnotes_math671_report.pdf}
\bibitem{PER001}
S.~Perlmutter, B.~P.~Schmidt, A.~G.~Riess, {\it The Accelerating Universe}, Scientific Background on the Nobel Prize in Physics 2011\\
{\scriptsize
\url{http://www.nobelprize.org/nobel_prizes/physics/laureates/2011/advanced-physicsprize2011.pdf}
}
\bibitem{SCH001}
A.~Schadschneider, {\it Einführung in die Hydrodynamik und ihre modernen Anwendungen}, Vorlesungsskript Wintersemester 2003/04, Universität Köln\\
\url{http://www.thp.uni-koeln.de/~as/Mypage/PSfiles/Hydro/hydro.pdf}
% Cooling table
% Sutherland & Dopita (1993, ApJS, 88, 253)
\end{thebibliography}
\end{document}
%%%%%%%%
% \subsection{Discretisation of the one dimensional Euler equation}
% \begin{eqnarray*}
% \partial_t\mathbf u(t,x)+\partial_x f(\mathbf u(t,x))=0
% \end{eqnarray*}
% \begin{eqnarray*}
% \mathbf u(t,x)=\begin{pmatrix}u_1(t,x)\\u_2(t,x)\\u_3(t,x)\end{pmatrix}=\begin{pmatrix}\rho(t,x)\\\rho(t,x)v(t,x)\\\epsilon(t,x)\end{pmatrix}
% \end{eqnarray*}
% \begin{eqnarray*}
% \mathbf w(t,x)=\begin{pmatrix}w_1(t,x)\\w_2(t,x)\\w_3(t,x)\end{pmatrix}=\begin{pmatrix}\rho(t,x)\\v(t,x)\\p(t,x)\end{pmatrix}
% \end{eqnarray*}
% \begin{eqnarray*}
% f(\mathbf u(t,x))=\begin{pmatrix}\rho(t,x)v(t,x)\\
% \rho(t,x) v^2(t,x)+p(t,x)\\
% (\epsilon(t,x)+p(t,x))v(t,x)
% \end{pmatrix}
% \end{eqnarray*}
%%%%%%%%