410-505 3.2 Hyperbolic Wave Equations Fall 2011

Home   |   Syllabus   |   Tools   |   Topic 1   |   Topic 2   |   Topic 3   |   Topic 4   |   Topic 5

This Page:   Home   >>   Topic 3   >>   Application 2

The Linear Wave Equation

The wave equation

where is the wave speed, and is a given source function, is a typical hyperbolic PDE.

This equation has a unique solution if

The one-dimensional wave equation

where is the wave speed, is the simplest example of a hyperbolic equation.

The wave equation on the infinite interval actually factorizes into simpler first order equations:

The equation is linear, and the solution can be expressed as a superposition of left-moving and right-moving waves

where the left-mover function and the right-mover function are determined by the initial conditions. These left- and right-movers are waveforms which propagate independently of one another with constant velocities , and without changing their respective shapes! On a finite interval the left- and right-movers cannot in general be chosen independently: they must cooperate to satisfy the boundary conditions. For example, Dirichlet boundary conditions can be satisfied by a standing wave, which has both left- and right-moving components.

The Simple Advection Equation

The right-mover equation

which has the analytical solution

where is the initial condition at , is an example of an advective equation.

Advection is a process in which a property at some point and time is determined by the shape of the function as it is carried past that point. This behavior is typical of fluid flow, for example: the density of the fluid at a point in space is determined by the density of the fluid element that happens to flow past that point at that time.

Contrast this with convection in which the density of the fluid element is determined by its position in space: hotter fluid is less dense and rises, and colder fluid is more dense and sinks; and the density changes as the fluid moves.

Flux-Conservative Equations in One Dimension

The simple advective equation above is an example of a flux-conservative equation

where is a vector of functions, and the vector is called the conserved flux of . These one-dimensional PDE's are discussed in Numerical Recipes Section 19.1.

Suppose that represents the density of fluid at point at time . The total amount (mass) of fluid between the boundaries is

The rate of change of fluid in the region is

which shows that can be interpreted as the fluid flux, that is, the mass of fluid flowing in the positive direction per unit time across an imaginary boundary located at and at time . The fluid moves such that its mass is conserved.

FTCS Algorithm for the Advective Equation

One possible discretization of the simple advective equation is obtained by introducing a spatial lattice

and a time step so that . Let .

The FTCS marching algorithm with this discretization is

Note that the spatial derivative has been approximated by a symmetric difference

Unfortunately, it turns out that this algorithm is unstable no matter how small the time step is made! This can be shown using von Neumann stability analysis. A mode

with wave number gets amplified by a factor , which is easily computed from the recurrence

Since for any mode, the FTCS algorithm is unconditionally unstable.

The Lax and Lax-Wendroff Methods for the Advection Equation

Lax modified the FTCS recurrence as follows:

It is easy to see that the mode amplification factor is

Note that if , which implies exact flux conservation for each mode! If modes will grow, making the solution unstable. However, if the modes will all damp to zero, which also does not conserve flux! So for the Lax method to give reasonable results, should be chosen close to and a little smaller so that large modes are damped.

The condition for stability

is the famous Courant-Friedrichs-Lewy criterion. The ratio on the left side of this condition is called the CFL number, and is especially important in computational fluid dynamics.

The criterion has a simple interpretation show in the figure. Consider the domain of dependency, shown shaded, which consists of all points in the past from which information can propagate at or slower than the wave speed . For any differencing scheme, the differencing domain consists of the set of points used to determine the next value of the solution.

The Lax-Wendroff method attempts to improve the solution accuracy by including terms of order in the discretization:

Here we have used .

The stability criterion for the Lax-Wendroff method is discussed in Numerical Recipes: it turns out to be the same Courant-Friedrichs-Lewy criterion as in the Lax method. The Lax-Wendroff method is considerably more accurate than the Lax method for the same time step .

Note that the added term of has the form of a diffusion term in the diffusion equation

which discretizes as

In general, a diffusive term in a recurrence formula has a damping effect on the solution amplitude.

Codes to Solve the Advection Equation

The instabilities of the FTCS, Lax, and Lax-Wendroff methods are demonstrated by running the codes:

Nonlinear Waves and Burgers Equation

J.M. Burgers, Adv. Appl. Mech. 1, 171 (1948), introduced the equation

as a simple model of shock propagation. This is basically a Navier-Stokes equation in one dimension without a pressure term. The convective term on the left is nonlinear. The diffusive term on the right represents the effects of viscosity.

The development of a shock can be seen by letting the kinematic vicosity . This gives the inviscid Burgers' equation

Compare this with the linear equation

where is a constant. The linear equation has the solution

where is any differentiable function. This solution represents a wave form with shape moving to the right with constant speed .

Now, in the inviscid Burgers' equation, the ``speed'' , i.e., the instantaneous speed of the wave form is proportional to its amplitude . This implies that a peak in the wave travels faster than a trough, which implies that the wave will tend to break. This is not allowed mathematically beacuse breaking implies that the solution becomes multiple valued. What actually happens is that a shock front develops: this is a moving point at which the solution is discontinuous.

The viscous term in Burgers' equation has two effects. First, it causes the wave amplitude to damp to zero in a diffusive fashion. Secondly, it prevents the development of a mathematical singularity at the shock front: the amplitude is continuous albeit varying very rapidly through the front.

Godunov's Method

This type of scheme was introduced by S.K. Godunov, Mat. Sb. 47, 271 (1959). This is an upwind differencing scheme which makes use of the solution to a local Riemann problem.

A Riemann problem is an initial value problem for a partial differential equation with a piecewise constant initial value function, i.e., the solution at is a step function. A Riemann solver is an exact or approximate algorithm for solving a Riemann problem.

The basic formula for updating is

where represents the average flux on the cells to the right and left of the lattice point respectively. These average flux values are computed from Riemann problems in the cells to the right and left of using upwind initial data

The solution to the Riemann problem on the left cell is

and for the cell on the right

Codes to Solve Burgers Equation

Euler's Equations of Hydrodynamics

Sod's Shock Tube Problem

The equations of fluid dynamics are mathematical statements of three fundamental physical principles:

The one-dimensional equations for the fluid dynamics of a gas can be written in conservation form as follows:

where is the density of the fluid, is the fluid velocity, is the energy per unit volume (length), and is the pressure. We need one more equation to close the system. This is the equation of state

where is the adiabatic gas index. For an ideal gas .

These equations can be written in vector form

where

The shock tube problem

A simple one dimensional model of a gas was introduced by G.A. Sod, J. Computational Physics 27, 1 (1978), to test the ability of various algorithms in solving fluid dynamics problems with shock wave behavior.

Sod considered a one-dimensional tube of unit length and the following initial conditions at :

This initial state can be produced by having a diaphragm in the middle of the tube. The gas to the left and right of the diaphragm is initially at rest. The pressure and density are discontinuous across the diaphragm. At , the diaphragm is broken. Two types of singularities then propagate through the gas:

To simulate a closed tube, reflection boundary conditions can be applied at . The shock tube then exhibits interesting behavior with shock waves and contact discontinuities bouncing back and forth in the tube and interacting with one another.

Godunov Methods and Riemann Solvers

Among the most interesting and difficult problems in computational fluid dynamics is the simulation of discontinuities like shock fronts. Simple finite difference schemes cannot handle this type of singular behavior.

Following the work of Godunov, Mat. Sb. 47, 271 (1959), which was based on his Ph.D. thesis, many effective shock-capturing schemes were developed for applications in astrophysics and the aerospace industry.

The linear advection equation

Consider the simple linear equation

where is a constant with dimensions of speed. Given an initial profile , the solution of this equation is easily seen to be , i.e., a waveform which moves at constant speed without changing its shape.

The Riemann problem

A simple form of initial condition is a step function or piece-wise constant value for , for example as shown in the figure. This type of initial condition defines a Riemann problem. Physically, this initial condition represents a shock front which moves with constant speed without changing its shape.

Even though this is such a simple problem with a simple solution, it is very difficult to simulate numerically. The reason for this is that the derivative is infinite at the discontinuity: mathematically it is a delta function. Most finite difference schemes assume that the solution is smooth, i.e., the derivatives are bounded, so that a Taylor series expansion in the spatial step size is valid. When this assumption is violated by a discontinuity, a first order scheme tends to smear out the discontinuity, and including higher orders results in unstable oscillations of the solution at the position of the discontinuity.

Integral form of the conservation law

To solve this problem, Godunov used the conservation form of the advection equation

where is the flux of the field .

The figure shows a few lattice sites on the space-time grid , that will be used to solve the problem numerically. If we consider the pair to be a vector function in the plane, then the conservation equation

is the divergence of the vector. Let us integrate this divergence over the rectangular region shown in the figure and use Gauss' integral formula to convert it to a line integral around the perimeter:

where the integrand in the line integral is the normal component of the vector field on the perimeter of the rectangle. Let's define the integral averages of on the top and bottom sides of the rectangle

and the time integral averages of the flux along the left and right sides of the rectangle

The line integral can be written

Now, if we interpret as the value of the solution at grid point at time step , then the value of the solution at grid point at the next time step is given by the formula

What remains is to specify the conserved half-step fluxes .

Godunov's upwind scheme for half-step fluxes

Godunov's suggestion for determining the half-step fluxes was to solve a pair of Riemann problems. For example, to determine , consider the Riemann problem on the interval for which is the center point

If the solution of this Riemann problem is denoted , then the Godunov flux is taken to be

For the linear advection equation, the solution of this Riemann problem is trivial

and hence

This gives an upwind scheme because if the waveform moves to the right and the left initial value covers the right boundary of the rectangular region; whereas if the waveform moves to the left and the right initial value covers the right boundary.

Substituting the Godunov flux values into the conservative update formula, we obtain the discrete solution

where the Courant-Friedrichs-Lewy number

This general Godunov approach can be applied to more complicated problems. For example, in Burgers' equation

and the current value of the solution, rather than the constant wave speed determines the choice of or in the equations above.

In the case of the 1-D Euler equations of gas dynamics, the solution and flux each have 3 components

Instead of a single wave speed, the solution to the Riemann problem involves finding the eigenvalues of a matrix: the solution involves several regions separated by left- and right-moving shock fronts and a contact discontinuity, instead of just a single shock front; the correct region must be chosen for to compute the Godunov flux.

Codes to solve the Euler Equations



© 2011   Richard J. Gonsalves

Department of Physics  |  University at Buffalo