Introduction

The method I chose for practice writing a higher dimension code is the 1st order accurate, split scheme, Godunov + Exact Riemann Solver to solve the 2D Euler equations for a cylindrical explosion. Given the higher dimension of the problem, new types of waves are present in the solution, namely shears. Shears are passively advected with the flow, as can be shown by combining the continuity equation and the corresponding momenta equations. This means that the velocity of the shear (i.e. velocity in the tangential direction relative to wave propagation) does not change across either shock or rarefaction. This means that the value of the tangential velocity only changes discontinuously across the middle wave of the Riemann problem, i.e. the contact discontinuity. The value of the shear velocity on the left side of the contact is just the tangential velocity of the left state, and vice versa. This makes for an extremely straightforward/easy modification to the exact Riemann solver.

Initialization

This method begins by initializing a 2D Cartesian grid, with a circle in the center. The primitive fluid variables (rho, x-velocity, y-velocity, pressure) are set inside and outside of this circle.

Here are some plots of the initial condition. You will see a circle of high density/pressure in the center of the mesh, with zero velocity in both x and y.

Density (pressure looks similar):

X-velocity (velocity in y is also zero):

The essence of split schemes

It can be proven for linear systems of equations of higher dimension, that splitting the problem into x and y components yields an exact solution. For non-linear systems, the effort leads to a highly accurate (but not exact) solution. The concept is as follows.

We begin with a 2D mesh of initial condition

\vec{u}^n = \{ \vec{u}^n_{i,j} \} ~\forall ~i,~\forall~j

where u is the vector of primitive variables. To solve the Euler equations for this given initial condition, we "split" the 2D Euler equations,

\frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho u \\ E \\ \rho v \end{pmatrix} + \frac{\partial }{\partial x} \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E+p) \\ \rho u v \end{pmatrix} + \frac{\partial }{\partial y} \begin{pmatrix} \rho v \\ \rho v^2 + p \\ v(E+p) \\ \rho u v \end{pmatrix} = 0

into 2, 1D problems. In the x-split direction, this is given by,

\frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho u \\ E \\ \rho v \end{pmatrix} + \frac{\partial }{\partial x} \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E+p) \\ \rho u v \end{pmatrix} = 0

and in the y- split direction this is,

\frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho v \\ E \\ \rho u \end{pmatrix} + \frac{\partial }{\partial y} \begin{pmatrix} \rho v \\ \rho v^2 + p \\ v(E+p) \\ \rho u v \end{pmatrix} = 0

Note that the first 3 equations of these different splits form just the 1D problem we have solved before. From a computational perspective, the difference for solving these sets of equations then is 1) the extra momentum equation, and 2) the presence of tangential velocity components that factor into E.

Now, the x- and y- split equations are solved on the grid independently. First the grid is sweeped in the x-direction, where it goes cell by cell on a given row and solves the initial value problem (IVP) given by the x-split equation above + the initial condition prescribed for the mesh. Once a row is is traversed, the next is traversed, until all of the grid has been solved in the x-split direction. This, by the way, is done the "normal" way— using the Godunov conservative update formula in conjunction with the exact Riemann solver.

The solution data for the x-split mesh is then stored in a 2D array, and used as the initial condition for the y-direction sweep. This sweep proceeds as above, this time moving cell by cell along a given column solving the IVP given by the y-split Euler equations above (and initial condition that is the SOLUTION to the previous x-split sweep). Once a column is traversed, the next is until the entire grid is sweeped. The solution to this last sweep is the solution to the full IVP (evolved to t+dt).

This process repeats for as many dt's needed to get to the final time of the simulation. Note, the dt is the SAME for both the x- and y-sweeps.

The method of choosing the best dt is described in Toro, chapter 16.

Upsides / Downsides

The good aspects of this choice of method (split, Godunov + Exact Riemann Solver) are: 1) accurate resolution of shear waves, 2) simple to construct.

The downside is: 1) not higher order (if wanted to add this - code could become quite complex to take into account the addition of shears into the TVD algorithms)

Download the code here

I followed Toro chapter 16 in writing this code. A previous blog post on coding nuances can be found ​here. The mathematica notebook I used for these plots is located at erica@crunch:~/Coding_Project/multidimensional. Here are downloadable versions of the ​code, ​module, and ​data file.

Results for cylindrical explosion

Here is a movie of density using extrapolating BC's: