wiki:u/ehansen/buildcode

Version 10 (modified by ehansen, 13 years ago) ( diff )

Building a 1D Hydro Code


Introduction

I've created this page to document the important aspects of building a hydro code. It contains general information that I have learned from reading Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction by E.F. Toro. In the near future, it will contain useful information from my own personal experience as I build my first hydro code.

The Euler Equations

These are the so called Euler Equations. They are also known as the fluid equations in conservative-law form. These are conceptually simple and make intuitive sense.

Name of Law Formula
Conservation of Mass
Conservation of Momentum
Conservation of Energy

Where is mass density, v is velocity (only in x-direction for 1D), p is pressure, and E is total energy per unit volume. E is further defined as where e is the specific internal energy. The specific internal energy depends on the equation of state. For an ideal gas where is the ratio of specific heats. Basically, these laws state that in a given volume, the change in a conserved quantity must be equal to the flux through the boundaries of that volume. In other words, the conserved quantity is in front of the time derivative and its flux is in front of the spatial derivative.

The Euler equations can also be written in integral form for a general control volume . To simplify the notation, let U be a column matrix containing the conserved quantities and F be the matrix of the corresponding fluxes. Now the Euler equations are:

In this form, it might be easier to see that the change in conserved quantity is equal to the flux coming into the volume minus the flux going out of the volume.

Solving the Riemann Problem

The Riemann problem is essentially the Euler equations with discrete initial conditions. Initial data has a left state and a right state, separated by a discontinuity.

Discretisation

Discretisation is the process of taking a domain and dividing it into cells. Each cell contains data corresponding to a specific problem initialization thus creating many Riemann problems. In other words, each neighboring pair of cells is now a "local" Riemann problem. Given two neighboring cells, the left cell is the left data state, the right cell is the right data state, and the intercell boundary is the discontinuity position. Now each local Riemann problem can be solved and that solution can be used to calculate the fluxes required to update the cells to the next time step.

Godunov's First Order Upwind Scheme

First, we want to use the integral form of the Euler equations (1). Next, we define a cell average for cell i at time step n + 1:

Where is the global solution as opposed to the local solution . However, we already know that the global solution can be written in terms of its fluxes. You just use equation (1) with a specific control volume:

The trick is to write the global solution in terms of the local solution. This can be done through a transformation of coordinates:

Where

Now we can combine (4) and (5) to get:

Substitute (6) and (7) into (3), and then that result into (2) gives:

Just looking at the right hand side of (8)…the first term can be rewritten just by using the definition of cell average that we defined in (2), and since the flux integrands are now constants, those integrals can be simplified. The final equation reads:


Program Outline

These are the basic steps in my program SEEQUOD (Solver for Euler EQUations in One Dimension). SEEQUOD currently uses an exact Riemann solver, but will later be updated to use approximate Riemann solvers as well. Also, SEEQUOD currently uses a first order Godunov upwind scheme, but again will later be updated to use higher order update schemes. These steps should be general enough to apply to many codes that are used to solve the Euler Equations.

  • Read all input data necessary for problem
  • Initialize domain cells
  • Begin time stepping procedure
    • Apply boundary conditions for "ghost" cells
    • Impose CFL condition to get appropriate time step
    • Begin flux calculation procedure
      • Solve local Riemann problem
      • Use solution from Riemann solver to compute fluxes
      • Repeat flux calculation procedure until fluxes are computed for all cells
    • Use fluxes to update conserved variables
    • Use updated conserved variables to update primitive variables
    • Repeat time stepping procedure until final time is reached
  • Output data to a file


Test Results

Five different tests were run in accordance with the tests in Ch. 6 of Toro's aforementioned book. All tests were run with a domain length of 1.0 and 100 computing cells. For the first five time steps a CFL of 0.18 was used to give the various waves time to develop. For the time steps following, a more natural CFL of 0.9 was used to obtain more accurate results. Here is a table of all the input data for each test:

Test Number Discontinuity Position Output Time
1 0.3 0.2 1.0 0.75 1.0 0.125 0.0 0.1
2 0.5 0.15 1.0 -2.0 0.4 1.0 2.0 0.4
3 0.5 0.012 1.0 0.0 1000.0 1.0 0.0 0.01
4 0.4 0.035 5.99924 19.5975 460.894 5.99242 -6.19633 46.0950
5 0.8 0.012 1.0 -19.59745 1000.0 1.0 -19.59745 0.01

Density, velocity, pressure, and internal energy were plotted for every test. The following images show the numerical results (points) alongside the exact solution (line). All plots appear to be the same as those presented by Toro. Therefore, the codes SEEQUOD and SEEQUODexact must be correct.

Test 1

Test 2

Test 3

Test 4

Test 5

Attachments (34)

Download all attachments as: .zip

Note: See TracWiki for help on using the wiki.