Topic4 Computational Fluid Dynamics

Document Sample
Topic4 Computational Fluid Dynamics Powered By Docstoc
					                           Topic 4: Computational Fluid Dynamics

                                     Lecture 4-1: Burgers Equation

                                        Monday, March 15, 2010


1 The equations of fluid dynamics                                                                      2

2 One Dimensional Burgers’ Equation                                                                   4

  2.1 Lax-Wendroff algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   6

  2.2 Godunov algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   7


1   The equations of fluid dynamics

The equations of fluid dynamics follow from conservation of mass and momentum.

Consider a volume V inside the fluid. The mass of fluid in this volume is given by

                                                    ρ dV ,

where ρ is the fluid density. The rate at which this mass decreases is determined by the rate at which fluid
leaves the volume
                                             ρ dV = − ρu · dS ,
where u is the fluid velocity and the integral on the right is taken over the surface of the volume with dS
being a surface element with direction along the outward normal. Using the divergence theorem

                                         ρu · dS =           · ρu dV ,

we obtain the continuity equation
                                              +      · ρu = 0 .
Next, consider conservation of momentum
                                               ρ      =F,


which is just Newton’s equation of motion for an element of fluid with unit mass. The total derivative on
the left has two contributions
                                          du ∂u
                                             =      + (u · )u ,
                                          dt    ∂t
the first term on the right represents the change in fluid velocity with t at a fixed point in space, and the
second advective term represents the change in fluid velocity due to motion of fluid from neighboring points
in space.

The force density F has three contributions:

  • external or body forces acting on the fluid, for example the force of gravity
                                                  Fgravity = ρg ,
    where g is the acceleration due to gravity
  • pressure forces due to neighboring fluid elements

                                          −      p dS = −        p dV ,

                                                 Fpressure = − p ,
    where p is the fluid pressure and the integrals are taken over the surface and volume of the element,
  • viscous forces due to internal friction or shearing stresses in the moving fluid
                                     Fviscous = µ       u + (µ + ξ) (   · u) ,

    where µ is the dynamic viscosity coefficient and ξ is the bulk viscosity coefficient of the fluid.

A special case that is interesting for many applications is that of incompressible flow
                                      ρ = constant ,         ·u=0.
Taking these forces into account results in the Navier-Stokes equations for incompressible viscous flow:
                                   ∂u                      1        2
                                      +u·      u=g−          p+ν        u,
                                   ∂t                      ρ
where ν = µ/ρ is the kinematic viscosity.

2   One Dimensional Burgers’ Equation

J.M. Burgers, Adv. Appl. Mech. 1, 171 (1948), introduced Burgers equation
                                            ∂u    ∂u   ∂ 2u
                                               +u    =ν 2 ,
                                            ∂t    ∂x   ∂x
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 ν = 0. This gives the inviscid
Burgers’ equation
                                         ∂u       ∂u
                                              +u     =0.
                                          ∂t      ∂x
Compare this with the linear equation
                                             ∂u     ∂u
                                                 +c      =0,
                                             ∂t     ∂x
where c is a constant. The linear equation has the solution
                                           u(x, t) = f (x − ct) ,
where f is any differentiable function. This solution represents a wave form with shape f (x) moving to the
right with constant speed c.

Now, in the inviscid Burgers’ equation, the “speed” c = u, i.e., the instantaneous speed of the wave form
is proportional to its amplitude u. This implies that a peak in the wave travels faster than a trough, which
implies that the wave will tend to break.

Breaking of two-dimensional surface waves is of course very familiar, see for example, Hokusai’s Great Wave
Off Kanagawa.

In one dimension, breaking is not allowed mathematically beacuse breaking implies that the solution u(x, t)
becomes multiple valued. What actually happens is that a shock front develops: this is a moving point at
which the solution is discontinuous.

2 ONE DIMENSIONAL BURGERS’ EQUATION                                                  2.1 Lax-Wendroff algorithm

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.

The code burgers.cpp solves Burger’s equation using various algorithms.

2.1   Lax-Wendroff algorithm

The Lax-Wendroff algorithm is constructed in two steps. First, the time and convective derivatives are
expressed in terms of a flow function F as follows:
                           ∂u    ∂u ∂u ∂F                            1
                              +u    =    +    ,            F (x, t) = u2(x, t) .
                           ∂t    ∂x   ∂t   ∂x                        2

This is the form of a conservation equation with F representing the current of the quantity u.

Second, a Taylor series expansion in the time step τ of all varables is made and terms up to and including
O(τ 2) are retained, e.g.,
                                                         ∂u τ 2 ∂ 2u
                              u(x, t + τ ) = u(x, t) + τ    +      2
                                                                     + O(τ 3) .
                                                         ∂t   2 ∂t

2 ONE DIMENSIONAL BURGERS’ EQUATION                                                          2.2 Godunov algorithm

The resulting algorithm can be expressed as a two-step formula:
                             1 n                τ
                    u∗ 1 =
                                 uj + un −
                                                     Fj+1 − Fjn +
                       2     2                 2h
                             ντ 1 n                           1
                                      uj+1 + un − 2un + un + un − 2un
                                               j−1       j                               ,
                             2h2 2                            2 j+2     j   j+1

                                    τ   ∗         ∗       ντ
                     j     = un −
                               j       Fj+ 1 − Fj− 1 + 2 un + un − 2un .
                                                               j+1  j−1   j
                                    h      2        2      h

2.2   Godunov algorithm

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.

Godunov’s algorithm is based on exact solutions to the Riemann Problem

2 ONE DIMENSIONAL BURGERS’ EQUATION                                                        2.2 Godunov algorithm

              u(x, 0)


A simple form of initial condition is a step function or piece-wise constant value for u(x, 0), 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 c 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 ∂u/∂x is infinite at the discontinuity: mathematically it is a delta

2 ONE DIMENSIONAL BURGERS’ EQUATION                                                            2.2 Godunov algorithm

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 h 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.

Godunov’s formula for updating u is
                                        τ                    ντ
                          un+1 = un −
                           j      j        Fj+ 1 − Fj− 1 + 2 [uj+1 + uj−1 − 2uj ] ,
                                        h      2        2    h
where Fj± 1    represents the average flux on the cells to the right and left of the lattice point j respectively.

These average flux values are computed from Riemann problems in the cells to the right and left of j using
upwind initial data
                         (+)    uj if uj > 0       (−)    uj if uj < 0
                       uj =                      uj =
                                0 otherwise               0 otherwise

The solution to the Riemann problem on the left cell is
                                                      1 (+) 2 1 (−) 2
                                      Fj− 1 = max      (u ) , (u )               ,
                                            2         2 j−1 2 j
and for the cell on the right
                                                    1 (+)    2       1 (−)   2
                                    Fj+ 1 = max       u          ,     u             .
                                        2           2 j              2 j+1

REFERENCES                                                                                  REFERENCES


[Recipes-C19-5] W.H. Press, S.A. Teukolsky, W. Vetterling, and B.P. Flannery, “Numerical Recipes in C”,
    Chapter 19 §5: Relaxation Methods for Boundary Value Problems,