410-505 | 3.2 Hyperbolic Wave Equations | Fall 2011 |
This Page: Home >> Topic 3 >> Application 2 |
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 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.
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.
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.
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.
The instabilities of the FTCS, Lax, and Lax-Wendroff methods are demonstrated by running the codes:
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.
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
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
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.
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.
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.
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.
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 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.